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ABSTRACT 

Studies of strong gravitational lensing in current and upcoming wide and deep pho- 
tometric surveys, and of stellar kinematics from (integral-field) spectroscopy at in- 
creasing redshifts, promise to provide valuable constraints on galaxy density profiles 
and shapes. However, both methods are affected by various selection and modelling 
biases, whch we aim to investigate in a consistent way. In this first paper in a series we 
develop a flexible but efficient pipeline to simulate lensing by realistic galaxy models. 
These galaxy models have separate stellar and dark matter components, each with 
a range of density profiles and shapes representative of early-type, central galaxies 
without significant contributions from other nearby galaxies. We use Fourier methods 
to calculate the lensing properties of galaxies with arbitrary surface density distribu- 
tions, and Monte Carlo methods to compute lensing statistics such as point-source 
lensing cross-sections. Incorporating a variety of magnification bias modes lets us ex- 
amine different survey limitations in image resolution and flux. We rigorously test 
the numerical methods for systematic errors and sensitivity to basic assumptions. We 
also determine the minimum number of viewing angles that must be sampled in order 
to recover accurate orientation-averaged lensing quantities. We find that for a range 
of non-isothermal stellar and dark matter density profiles typical of elliptical galax- 
ies, the combined density profile and corresponding lensing properties are surprisingly 
close to isothermal around the Einstein radius. The converse implication is that con- 
straints from strong lensing and/or stellar kinematics, which are indeed consistent 
with isothermal models near the Einstein radius, cannot trivially be extrapolated to 
smaller and larger radii. 

Key words: gravitational lensing - stellar dynamics - galaxies: photometry - galax- 
ies: kinematics and dynamics - galaxies: structure - methods: numerical 



1 MOTIVATION 

1.1 Learning from galaxy density profiles and 
shapes 

Observational constraints on galaxy density profiles and 
shapes can be used to study a great variety of problems, 
from basic cosmology to the connections between dark mat- 
ter and baryons. For example, the spherically-averaged form 
of the dark matter halo density profile, and the distri- 
bution of halo shapes, both depend on cosmological pa- 
rameters to the extent that those p arameters affect the 
process of structur e formation (e.g., iBuUock et al.l l200ll : 
lAllgood et aLll2006l ). and on the physics of galaxy forma- 
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tion to the extent that baryons modify the dark matter dis- 

I — — 1 1 — — ^ ' 1 ' ' 

tribution (e.g., Blumenthal et al. 1986; Gncdin ct al. 200: 



[Kazantzidis et al. 2004; Naab et al. 2007; Rudd ot al. 200S 
While there is still some disagreement in the literature about 
the theoretical predictions, even among those papers cited 
above, one key result is that there is considerable scatter in 
both the density profiles and shapes of dark matter halos, 
which presumably re flects different formation histories (e.g., 
IWechsler et"aLll2002l ). 

Several observational tools have been used to constrain 
the density profiles and shapes of galaxies, at different length 
and mass scales, and redshifts. The main tools we consider 
are gravitational lensing and stellar kinematics. X-ray data 
are also useful for understanding the density profiles and 
shapes of massive galaxies and galaxy clusters, but are be- 
yond the scope of our investigation. 

Gravitational lensing is the deflection of light from dis- 
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tant sources by the gravitational fields of intervening lens 
galaxies. Strong lensing occurs in the central regions of 
galaxies (and clusters) where the light bending is extreme 
enough to produce multiple images of a background source. 
It can be used to study the mass distributions of individual 
galaxies on projected scales of typically several kpc (for a 
review, see iKochaneklliooel ). Weak lensing is a complemen- 
tary phenomenon in which the deflection of light slightly 
distorts the shapes of background sources without creating 
multiple images (for a review, see iBartelmann fc Schneider] 
I2OOII ). Weak lensing probes galaxy mass distributions on 
scales from tens of kpc out to about 10 Mpc, but only yields 
constraints on ensemble averages, since currently weak lens- 
ing by galaxies can only be detected by stacking many lens 
galaxies. While galaxy clusters can by studied on an indi- 
vidual basis using weak lensing, our focus is on galaxies. 

High-quality (two-dimensional) data on the kinemat- 
ics of stars as well as gas in the inner parts (a few to 
ten kpc) of galaxies are now readily available, thanks in 
particular to strong progress in integral-field spectroscopy 
in the last decade. Kinematics in the outer parts of late- 
type galaxies can often be observed from the presence 
of neutral hydrogen (e. g., Bosmal 198ll: van Albada et al.l 
1 19851 : iPersic et al. I Il996l : iNoordermeer et all |2007| ). In the 
outer parts of ear ly-type galaxies, however, cold gas is 
scarce (but see e.g. Franx et al. I I1994I : iMorganti et al.l1l997l : 
IWeiimans er aT 2008) , so we are left with discrete kinematic 



tracers sucn as planotary nebulae and g lobular cluster (e.g. 
ICote et al.ll2003l : iDouglas fc et al.ll2007r ) out to tens of kpc 
Current and upcoming wide and deep photomet- 
will reveal a vast nu mber of strong l ensin. 



ric surveys will reveal a vast number ot strong lensing 
events (e.g., [ Fassnach t et al.l 2004: Koop mans et al] |2004| : 
iKuhlen et al] 12004 : ,Marshall et al. ,2005 ) , and will also al- 
low for extensive, complem entary weak lensing analysis (e.g. 
iTvsonI [200^ : iKaised [2001 ). At the same time, there is a 
rapid increase in the availability of two-dimensional kine- 
matics of galaxies near by (e.g., lEmsellem fc et al] l2004l : 
McDermid fc et al.l200"^), and even at high(er) redshift (e 



van der Marel fc van Dokkum|[2007l : iBouche fc et al.|[200' 



These data sets, individually as well as combined, may 
provide strong constraints on galaxy density profiles and 
shapes, but the analyses involved are subject to selection 
and modelling biases. 

1.2 Selection and modelling biases 

If we want to use strong lensing and kinematics to make 
robust tests of cosmological predictions, we need to answer 
two questions: Can strong lensing and kinematic analyses 
yield accurate constraints on galaxy density profiles? Are 
the galaxies in which we can make the measurements (es- 
pecially in the case of strong lensing) representative of all 
galaxies? We refer to these two concerns as modelling biases 
and selection biases, respectively. 

The question of selection bias is particularly important 
given the diversity in the galaxy population. It is well known 
that strong lensing favours early-type over lat e-type galax- 
ies, and massive galaxies over dwarfs (e.g., [Turner et al.l 
1 19841 : iFukugit a fc Turneijll99il ). But even within the pop- 
ulation of massive early-type galaxies, to what extent does 
strong lensing favour galaxies whose dark matter halos have 
inner slopes that are steeper than average, or concentra- 



tions that are higher than average? Also, to what extent 
does strong lensing favour galaxies with particular shapes 
and/or orientations with respect to the line-of-sight? To 
phrase these questions formally, consider some parameter 
X describing the galaxy density profile or shape (e.g., the 
inner slope of the dark matter density profile). We need to 
understand how the distribution psL(a;) among strong lens 
galaxies compares to the underlying distribution p{x) for 
all galaxies. Any analysis that includes strong lensing with 
some other technique used to study the same systems will 
suffer from strong lensing-related selection biases, since we 
never get to choose which systems will be a strong lens. 
Selection biases are critical when we want to interpret con- 
straints on lens galaxy density profiles and/or shapes from 
strong lensing, possibly in combination with kine matics (e.g. 
iRusin fc Kochane^l2005l : iKoopmans et al.ll2006l ) in compar- 
ison with predictions from cosmological models. 

On the other hand, modelling biases may occur because 
of (often unavoidable) assumptions made when analysing 
gravitational lensing and/or kinematic data because of the 
finite number of constraints available from the data. Since 
galaxies are in general non-spherical, we can only correctly 
interpret the observations if we know the viewing direction, 
and even the n the deprojection might not be unique (e.g. 
iRvbickil Il987l ). Moreover, strong lensing is subject to the 
so-called mass-sheet degeneracy: part of the deflection and 
magnification of the light from the background source can 
be due to mass along the line-of-sight that is not associated 
with the lens galaxy itselfl- In addition, the constraints from 
strong lensing on the galaxy density are typically limited in 
radius to around the Einstein radius, and even then, without 
secure, unaffected measurements of the flux of the lensing 
images, the density profile is difficult to recover. Kinematics 
can be obtained over a larger radial extent, but in particular 
stellar kinematics suffer from the so-called mass-anisotropy 
degeneracy: a change in the measured line-of-sight velocity 
dispersion can be due to a change in total mass, but may 
also be the result of velocity anisotropy. 

In a series of papers, we are developing a pipeline for 
using realistic galaxy models to simulate strong lensing and 
kinematic data, and assess how both selection and modelling 
biases affect typical strong lensing and kinematic analyses. 
In this first paper, we present the simulation pipeline for 
point-source lensing. We focus on strong lensing of quasars 
by early-type, central galaxies at different mass scales, using 
two-component mass profiles (dark matter plus stellar com- 
ponent) that are consistent with exist ing photometry and 
stack e d weak lensing data from S PSS (iKauffmann fc et al.l 



I2OO3I : iMandelbaum et all l2006bl). A^-body and hydrody- 



dy, and 

namic simulations. In Mandclba um et alT (|2008l . hereafter 
Paper II), we use this pipeline to study selection biases in 
strong lensing surveys. In future work we will address mod- 
elling biases in strong lensing and kinematics (both when 
studied separately and when combined). Our overall goal is 
to determine how strong lensing, kinematics, and possibly 
other probes (weak lensing, and X-ray data for cluster mass 



^ Part of this mass-sheet degeneracy may be overcome in future 
surveys through the measurements of time delays for an adopted 
Hubble constant. 
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scales) can be used to obtain robust, unbiased constraints 
on galaxy density profiles and shapes. 

To make our presentation coherent, and to clarify our 
notation and terminology, we first review the analytic treat- 
ment of ellipsoidal mass distributions with different profiles 
and shapes (Section [5}, and the basic theory of strong lens- 
ing (Section[3]). We then describe our simulation pipeline for 
point-source lensing in depth. In Section |4] we present our 
choices for the masses, profiles, and shapes of the galaxy 
models. In Section [5] we discuss the numerical methods we 
use for lensing calculations, including numerous tests. A 
summary of the simulation pipeline, and some implications 
for strong lensing analyses are given in Section [6] 



2 GALAXY DENSITY PROFILES AND 
SHAPES 

In this section we describe our treatment of galaxy density 
profiles and shapes. While our approach is fairly conven- 
tional, it is important to present it carefully to clarify our 
notation and terminology, collect useful technical results, 
and provide a firm foundation for the work in this series of 
papers. This section focuses on the formal framework; the 
specific galaxy models used in our simulations are discussed 
in Section |4l 



2.1 Notation 

We compute distances using a fiat ACDM cosmology with 
f2m = 0.27 and h = 0.72. We use x, y, and z to denote in- 
trinsic, three-dimensional (3d) coordinates, and x' and y' for 
the projected, two-dimensional (2d) coordinates. Similarly, 
r is the intrinsic radius (r^ — x'^ + y^ + z'^) and R' is the pro- 
jected radius {R'^ — x'"^ -f y'^). For non-spherical galaxies, 
we use a, b, and c for the major, intermediate, and minor 
semi-axis lengths of the intrinsic, 3d density profiles; and 
we use a' and b' for the major and minor semi-axis lengths 
of the projected, 2d surface densities. Subscripts "dm" and 
"★" are used to indicate whether a quantity describes the 
dark matter or stellar component of the galaxy model. Any 
gas component that has not cooled to form stars — which 
is expected to be subdominant rel ative to the stellar com- 
ponent for early type g alaxies (e.g. iRead fc Trenthamll2005l : 
iMorganti &: et aL 12006 ) — is implicitly included in the com- 
ponent that we label dark matter. 



2.2 Ellipsoidal shapes 

We consider mass density distributions p(x,y,z) = pirn) 
that are constant on ellipsoids 



2 

0^ c-' 



(1) 



with a ^ b c. The major semi-axis length a is a scale 
parameter, whereas the intermediate-over-major (6/a) and 
minor-over-major (c/a) axis ratios determine the shape. In 
the oblate or prolate axisymmetric limit we have a — b > c 
(pancake-shaped) or a > b = c (cigar-shaped), respectively, 
while in the spherical limit a — b — c (so then m — r/a). 
Note that m is a dimensionless ellipsoidal radius. 

Under the thin-lens approximation, the gravitational 



lensing properties are fully characterised by the mass den- 
sity projected along the line-of-sight. We introduce a new 
Cartesian coordinate system (x" ,y" , z"), with x" and y" 
in the plane of the sky and the 2"-axis along the line-of- 
sight. Choosing the a:"-a xis in the (x, y)-pla.ne of t he intrin- 
sic coordinate system (cf. Ide Zeeuw fc FraroJllOSgl and their 
Fig. 2), the transformation between both coordinate systems 
is known once two viewing angles, the polar angle i9 and az- 
imuthal angle ip, are specified. The intrinsic z-axis projects 
onto the y"-a.xis; for an axisymmetric galaxy model the y"- 
axis aligns with the short axis of the projected mass densityQ 
but for a triaxial galaxy model the j/"-axis is misaligned by 
an an gle ip € [—tt/2, tt/2] such that (cf. equation B9 of .Franxl 
I1988D 



tan 2-0 



T sin 2(/9 cos 
sin^ ■& + T (sin^ p cos^ ■& — cos^ (p^ 



(2) 



where T is the triaxiality parameter defined as T = (a^ — 
b'^)/{a? — (?)■ A rotation through ip transforms the coordi- 
nate system (x" ,y" , z") to {x',y',z') such that the x' and 
y' axes are aligned with the major and minor axes of the 
projected mass density (respectively) , while z' — z" is along 
the line-of-sight (see also Section [4.4.31 below) . 

Projecting p(m) along the line-of-sight yields a surface 
mass density T.{x',y') — E(m') that is constant on ellipses 
in the sky-plane. 



E(m') 



, ^ , / abc 
p(^)d. =_2 



p(m) m du, 



(3) 



where we have used z' — abc sinh(u)/(a'6') + constant, and 
m = m' cosh(it). The sky-plane ellipse is given by 



,2 x'^ ^ y'^ 



(4) 



The projected major and minor semi-axis lengths, a' and 
depend on the intrinsic semi-axis lengths a, b, and c and the 
viewing angles -d and as follows: 



2 A' 



2A^ 



B - VB^ - 4^2 ' B + VB'^ - 4^2 ' 

where A and B are defined as 



(5) 



= a^b^ cos^ i9 -|- (a^ sin^ ip + b'^ cos^ p)<? sin^ (6) 

B 



a (cos ip cos 1? -I- sin p) 
+b'^ (sin^ ip cos^ i? + cos^ p) + (? sin^ i9. 



(7) 



It follows that A = a'fe' is proportional to the area of the 
ellipse. 

The fiattening fe'/a' of the projected ellipses actually de- 
pends on the viewing angles (i9, p) and the intrinsic axis ra- 
tios (6/a, c/a), and is independent of the scale length. Since 
the intrinsic and projected semi-major axis lengths are di- 
rectly related via the left equation of (O, the scale length 
can be set by choosing either a or a'. 

Finally, we define the negative logarithmic slope of the 
mass density 7(m) = — din p(m)/dln m and of the surface 



2 For an oblate galaxy, the alignment = 0) follows directly 
from equation l[2]l since T = 0. For a prolate galaxy it is most 
easily seen by exchanging a and c (c > 6 = a), so that the z-axis 
is again the symmetry axis (instead of the a;-axis). 
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mass density 7'(m') — — d In E(m')/d In m'. Substituting 
equation ([3}, the latter follows as 

'/ 'N f^['yim)-l]pim)m' du 

-f{m) = -^ , (8) 

Jg p^m) m au 

which in general has to be evaluated numerically. 



2.3 Choice of density profiles 

We consider two families of density profiles, motivated by 
both observations and simulations of galaxies. Historically, 
there has been considerable interest in cusped density pro- 
files that have a shallower power-law at small radii and a 
steeper power-law at large radii, wit h a smoot h tran sition 
in between. This family includes the iHernauisd (|l990[ ) pro- 
file that is commonly used to model the stellar components 
of early-type galaxi es and spiral galaxy bulges, along with 
(generalised) NFW (Navarro et al.j |l997l ) profiles often used 
to describe the dark matter profiles of simulated galaxies. 

An alternative, observationally m otivate d fam ily of 
models is obtained by deprojecting the ISersid (|l968l ) pro- 
file that is known to give a good fit to the surface bright- 
ness profiles of early-type galaxies and spiral galaxy bulges. 
There have been claims that deprojected Sersic profiles ac- 
tually provide a better fit to simulated dark matter ha- 
los than (generalised) NF W profiles (e.g., iNavarro fc et al] 
l2004l : iMerritt et ai]|2005l ) . Several recent studies argue that 
a Ser sic profil e that is not deprojected — also referred to 
as an [ Einasto (Il965l) profile — provides a n even better fit 
(|Merritt et al. |2006| ; iGaol IL. and Navarrol ) , although when 
compared against a generalised NFW or Sersic profile the 
improvement is margi nal at best. 

X -ray (e.g., .Allen et al.l l200ll . [20o3: IVikhlinin et all 

20051. l2006t) and we ak lensing (e.g.^ iKneib et al.l l2003l : 



Limousin et al.ir2007l ) observations of individual galaxy clus- 
ters suggest that outside of ~ 50 kpc, the density profiles 
are consistent with the NFW profiles seen in N-body simu- 
lations (no attempts were made to compare with the Sersic 
model, however). The mass range for these observations is 
about a factor two above our higher mass model. 

Since there is still no consensus opinion about which 
models are best, we consider how the choice of density profile 
affects our conclusions by using both cusped density profiles 
and deprojected Sersic density profiles 

2.4 Cusped density profiles 

The cusped density distribution with inner slope 7 and outer 
slope n, 

Po 



p(m) = 



(9) 



mT(l + m)"-T ' 

includes the following well-known spherical (m = r/a) den- 
sity p rofiles: (7, n) = (1, 4), the Hernquist p rofile (i fernquistl 
Il99d ): (2,4), the Jaffe profile jjaffelfigsl ): and (1,3), the 
NFW profile (iNavarro et al.lll99it ). 

The mass enclosed within the ellipsoidal radius m (for 
inner slopes 7 < 3) is 



M(m) 



1 + m 



3-7 



X 2F1 



n — 7, 3 — 7; 4 — 7 



'1+m 



(10) 



where 2^1 [a, /3, 7; a;] is the hypergeometric function. For 
outer slopes n > 3, equation (|10|) reduces to Mijn) = 
Anabcpo /3[n — 3, 3 — 7; m/ {l-\-m)]. When m — > 00, the latter 
incomplete beta function becomes the complete beta func- 
tion /3[n — 3, 3 — 7], and we obtain a finite total mass. 

Substituting n = 4 into equation [TO] yields M (m) = 
47ra6cpo['T^/(l+'T^)]'^~^/(3— 7), so that for 7 = 1, we find the 
total mass M — 2nabcpo of the Hernquist profile, which we 
adopt for the stellar component when using cusped density 
profiles. For outer slopes n ^ 3, the total mass is infinite, 
but for certain half-integer values of 7 and n the expression 
for the enclosed mass simplifies significantly. For example, 
for n = 3, the value which we adopt for the outer slope of 
the dark matter profile, we find 



M(m) 
AzTTabcpo 



m(2 + 3m)/[2(l + m)2], 

- 2V^(4m + 3)/[3(l + mf^% 
m/(l -I- m). 



ln(l + m 
2sinh"^( 
= < ln(l +ra) 

2 sinh"^ ( Vm) - 2^m/(l+m) 
\n{l + m), 



(11) 



for inner slope values of 7 = {0, 1/2, 1, 3/2, 2}, respectively. 

Although for certain values of 7 and n, lengthy analytic 
expressions for the surface mass density E(m') can be de- 
rived, we evaluate the integral in equation ((3} numerically. 
The logarithmic slope 7(m) of the cusped density profile is 



7(m) = (7 + mn)/{l + m), 



(12) 



while j'{m') of the surface mass density follows from rm- 
merical evaluation of equation (|8]) . 



2.5 Deprojected Sersic density profiles 

It has long been known that the surface brightness profiles 
of early - type g alaxie s and of spiral galaxy bulges are well 
fit by a lsirsid |l968") profile 7(7?) oc exp[-(7i/7?e)^^"], with 
the effective radius Re enclosing half of the total light. A 
key conceptual difference from the cusped models is that 
the profile does not converge to a particular inner slope on 
small scales. The deprojection of the Sersic profile has to be 
done numerically. Howev er, the analytic density profile of 
IPrugniel fc SimienI (|l997l ). 



pirn) = exp m^''"j , 



(13) 



provides a good match to the deprojected Sersic profile when 
the inner negative slope is given by 



, 0.6097 0.05563 

Pn = 1 \ ^ . 



(14) 



n 

The enclosed mass for the Prugniel-Simien model is 

M(m) = 4Tvabcpo nbl^"--^^" ^[{3 - p„)n,b„m^^"], (15) 

where 7[p;a;] is the incomplete gamma function, which in 
the case of the total mass reduces to the complete gamma 
function r[p] — 7[p; 00]. 

The expression for the surface mass density is, to high 
accuracy, the Sersic profile 



E(m') = Eo exp j^— 6„ (m')^''"j 



(16) 
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Given the enclosed projected mass 

M' (rn) — 27ra'6'Eo nb^^" 'y[2n,bn {m')^^" 



(17) 



the requirement that the total intrinsic and projected mass 
have to be equal yields a normalisation 



En = 



abc 2r[(3-p„)n] 



2n\ 



(18) 



The value of &„ depends on the index n and the choice 
for the scale length. The latter is commonly chosen to be 
the effective radius Re in the surface brightness profile, 
which contains half of the total light. We adopt a similar 
convention requiring that the ellipse m' — 1 contains half 
of the projected mass. This choice results in the relation 
r[2n] = 2 7[2n, which to hig h precision can be approxi- 
mated by (|Ciotti fc Bertinll 19991 ) 



14 1 46 1 

6,1 = 2 n \ 1 . 

3 405 n 25515 



(19) 



The logarithmic slope of the Prugniel-Simien and Sersic 
profile are related by 

jim)=p„ + 'y'{m), (20) 
where the logarithmic slope of the surface mass density is 
7'(m') = (fen/n)(m')'/". (21) 



3 STRONG LENSING 

In this section, we review the stron g lensing concep t s that 
are most impor t ant fo r our work. See fSchneider et al.l (Il992l ) 
and iKochanekl lj2006h for more discussion of strong lensing 
theory. The specific lensing calculations used in our simula- 
tions are discussed in Section [5] 



3.1 Basic theory 

The gravitational lensing properties of a galaxy with sur- 
face mass density E(i?') are characterised by the lens po- 
tential <j) (in units of length squared) that satisfies the two- 
dimensional Poisson equation 



V (j) = 2k , where k = E/Ec 



(22) 



Here « is the surface mass density scaled by the critical 
density for lensing (Ec, see Section [3.2p . and is referred to 
as the "convergence." The positions and magnifications of 
lensed images depend on the first and second derivatives of 
the lens potential, respectively. The image positions are the 
solutions of the lens equation, 



(23) 



where a = \/<j} is the deflection angle (in units of length), 
while 6 = {61,92) and /3 = (/3i,/32) are two-dimensional an- 
gular positions on the sky in the lens and source planes, 
respectively. The angular coordinates are related to the pro- 
jected physical coordinates by 0i = x' /Dl and 62 = y' /Dl, 
where Dl is the angular diameter distance to the lens; and 



/9i = x' /Ds and (32 = y' /Ds, where Ds is the angular diam- 
eter distance to the sourcelf] The magnification of an image 
at position {x',y') is given by 



dx'2 



1 - 



9^ 



dx'dy' 



(24) 



A typical lens galaxy has two "critical curves" along which 
the magnification is formally infinite, which map to "caus- 
tics" in the source plane. The caustics bound regions with 
different numbers of images (see Section [5.71 for examples). 

If the lens is circularly symmetric, the deflection is ra- 
dial and has amplitude 



a{R') 



2 

R' 



R' 



k{R') R' dR' . 



(25) 



The outer or "tangential" critical curve corresponds to the 
ring image that would be produced by a source at the ori- 
gin. The radius of this ring is the Einstein radius -Rein, 
which is given mathematically by the solution of the equa- 
tion a(-Roin) = Rein- This definition is equivalent to saying 
that the Einstein radius bounds the region within which the 
average convergence is unity, such that the enclosed mass 



TT-Rein^c By contrast, the inner or "radial" critical 



curve has radius i?rad given by the solution of da/dR' — 1. 

If the lens is non-circular, the conditions for the critical 
curves and caustics are more complicated, but we can still 
retain some concepts from the circular case. In particular, we 
can say that the outer/tangential critical curve is principally 
determined by the enclosed mass, whereas the inner/radial 
critical curve is determined by the slope of the deflection 
curve. When we map the critical curves in the lens plane to 
the caustics in the source plane, the arrangement of curves 
is inverted: the outer critical curve corresponds to the inner 
caustic; while the inner critical curve corresponds to the 
outer caustic|3 Since the outer caustic bounds the multiply- 
imaged region, it determines the total lensing cross-section. 
This means the cross-section is sensitive to the surface mass 
density (in particular its logarithmic slope) at small radii 
in the lens plane, which will be important to bear in mind 
when interpreting our results. 

While a non-circular lens has non-circular critical 
curves, it is still occasionally useful to characterise the scale 
for strong lensing with a single Einstein radius. We gener- 
alise the deflnition of the Einstein radius to a non-circular 
lens by using the monopole deflection, or the deflection an- 
gle produced by the monopole moment of the lens galaxy, 

1 /■27r rR' 

ao(R')^—^ / k(R,§) RdRde ^ R'k(R'), (26) 
T^R Jo Jo 

where k{R') is the average convergence within radius R' , 

Note that because of this trivial conversion between angular 
and physical units, we do not explicitly distinguish between angles 
and lengths in the lens plane. For example in the text we refer 
to the deflection angle in physical units (e.g., by stating that 
Q(Kcin) = ^ein) whereas our plots may show a in angular units. 
* If the lens is highly flattened, the tangential caustic can actually 
pierce the radial caustic; some specific examples are discussed 
below (e.g., Section l5.7| l. However, this situation is unusual, and 
it does not substantially alter the present discussion. 
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and {R, 9) are polar coordinates on the sky. The Einstein 
radius is then given by ao{Rcin) ~ Rein, or equivalently 
K{Rein) = 1- This is the definition that emerges naturall y 
from models of non-spherical lenses (e.g.. lCohn et al.]|200l[ ). 
and it matches the conventional definition in the circular 
case. 

It is instructive to think of lensing in terms of Fermat's 
principle and say that lensed imag es form at stationary 
point s of the arrival time surface ('e.g.. lBlandford fc NaravanI 
1 19861 ). Since the arrival time is a 2d function of angles on 
the sky, there are three types of stationary points: minima, 
maxima, and saddle points. We can classify them using the 
eigenvalues of the inverse magnification tensor. 



A± = (1 — k) ± 7 , where 7 = [(1 — k)^ — jj. 



-ni/2 



(27) 



Here, 7 is the sheai[f|, which quantifies how much a re- 
solved image is distorted. The eigenvalues are both positive 
for a minimum, both negative for a maximum, and mixed 
for a saddle point. Since our galaxies have central surface 
densities shallower than S oc R'~^, the lensing "odd im- 
age theorem" applies: the total number of images must be 
odd, and the number of minimum, maximum, and saddle 
point images must sat isfy the following relation (|Burkell98ll : 
[Schneider eralll 19921 ) : 



iVsad + 1 . 



(28) 



In practice, images at maxima are rarely observed because 
they form near the centres of lens galaxies and are highly 
demagnified; the only secure observation of a maximum im- 
age in a galaxy-scale len s required a deep and dedicated 
search l|Winn et al.ll2004 ). We do compute maxima because 
they are useful in checking our numerical methods (see Sec- 
tion [5]2]), but we focus our analysis on minimum and saddle 
point images because they constitute the bulk of lensing as- 
trophysics. 



3.2 Redshift dependence 

The redshifts of the lens [zl) and source [zs) play an impor- 
tant role in determining what region of a galaxy is relevant 
for strong lensing. Roughly speaking, strong lensing occurs 
in the region where the surface mass density exceeds the 
critical density for lensing, 

Ds 



- Ai,GDlDls ' ^^^^ 
where the D values are the angular diameter distances to 
the lens (-Dl), to the source {Ds), and from the lens to the 
source (Dls)- Note that Dls can be computed for a flat 
universe using the relation 



Dls = Ds 



l + ZL 



Dl 



(30) 



1 + zs, 

and for more general cosmologies, see lHog3 l|l999l ). For a 
fixed source redshift, Ec is lowest when the lens is roughly 



^ In general the shear has two components conve niently expressed 
in complex notation (e.g. [Schneider et al.lll993) , but as usual in 
strong lensing analyses we only use the amplitude (or norm) for 
which 7 is the standard symbol. The context should make clear 
whether we are using 7 to refer to a lensing shear or to the neg- 
ative logarithmic slope of the density profile. 



halfway between the observer and source, and it increases 
as the lens moves toward the observer or source. For a fixed 
lens redshift, the critical density is formally infinite for zs ^ 
Zl, and it decreases monotonically as the source redshift 
increases past zl- 

We adopt fiducial redshifts of zl = 0.3 for the lens 
and Zs = 2 for the source. The resulting critical density, 
Ec = 2389M0pc~^, is such that quasar images generally 
appear at about one effective radius from the centre of the 
lens galaxies. At this radius, the stellar component and dark 
matter halo may both play significant roles in the lensing 
signal, which is important to keep in mind when interpreting 
our results in this and subsequent papers. 

Our results apply equally well to any other combination 
of lens and source redshifts that yield Ec — 2389 M© pc~^. 
The reason is that if we fix the galaxy density profile and 
vary zl and zs so as to keep Ec fixed, the strong lensing 
region will always have the same physical scale and hence 
the same relation to the density profile. If we instead vary 
Zl and zs in a way that increases Ec, that will tend to push 
the lensed images to smaller radii and hence make the stel- 
lar component somewhat more dominant. Conversely, vary- 
ing the redshifts in a way that lowers Ec will tend to make 
the dark matter component more significant. We examine 
the effect of varying zl and zs on our results explicitly in 
subsequent papers. 

We note that the lens redshift also affects strong lensing 
in the conversion from physical to angular scales. Even if the 
physical size of the lensing region stays fixed as we vary the 
source and lens redshifts, the angular separation between 
the lensed images will vary with zl- This may lead to 06- 
servational selection effects: for example, resolution-limited 
surveys may miss lenses for which zi, is high and the images 
cannot be resolved, while spectroscopic surveys may miss 
lenses for which 21, is low and the images fall outside the 
spectroscopic slit or fiber. However, our focus in this work is 
on physical selection biases in strong lensing, which means 
that we consider the lens and source redshifts only in terms 
of how they affect Ec. 

Our choice zs ~ 2 for the source redshift is typical 
of quasar lenses in the CASTLE^ sample, but our choice 
Zl ~ 0.3 places the lens galaxy on average at a smaller dis- 
tance. The latter is, however, closer to the typical lens red- 
shift in the SLAGS sampl e of galaxy /galaxy strong lenses 
("Bol ton et"all |2006| . l2008al ). but those lenses are exphcitly 
selected to have source redshifts zs < 0.8 (|Bolton et al.l 
|2004 ). As a result, our fiducial value of Ec is typically be- 
tween that of CASTLES and SLAGS lenses. The Einstein 
radii at ~ 1 J?c in our c ase appear at larger r adii in the CAS- 
TLES sample (2- 3 fi>.: lRusin fc et al1l2003h. and at smaller 



radii in the SLACS sample (0.3 — 0.9 R^; iKoopmans et all 
I2OO6I ). The statistics of the SLACS survey are much more 
involved than the statistics of quasar lens surveys, because 
SLACS lenses have extended rather than point-like sources, 
and the survey relies on fiber spectra from S PSS. SLACS 
lens statistics have recently been examined bv lDobler et al.l 



« CASTLES (see |http://cfa^ www.harvard.edu/castles7) is a col- 
lection of uniform HST observations of mostly point-source lenses 
from several samples with differing selection criteria, rather than 
a single, uniformly-selected survey. 
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(|2008l l. and we do not consider them explicitly in the current 
work. 



4 CONSTRUCTION OF GALAXY MODELS 

In this section we construct the galaxy models used in our 
simulations. We seek models with a variety of density pro- 
files and shapes that not only are realistic but also sample 
the range of systematic effects in strong lensing. To that 
end, we must carefully consider the mass and length scales 
(Sections 14. 1144. 2|) . the profiles (Section 14. 3p . and the shapes 
(Section 14.41) of both the stellar and dark matter compo- 
nents. In addition to a main set of simulations with cusped 
density profiles, we also create a subset of simulations based 
on deprojected Sersic density profiles (Section [43} . After 
discussing the choice of the model parameters, we briefly de- 
scribe the part of our pipeline that computes surface mass 
density maps for these models (Section l4.6|l . which can then 
be used for lensing calculations (Section [S]). 



4.1 Mass scales 

We use stellar and dark matter halo masses from SDSS for 
early-type galaxies, where analyses of the spectra gi ve an 
estimate of the stellar mass (Ka uffmann fc et al.|[2003l '). and 
weak lensing analyses from 0.02-2fe~^ M pc yield the dark 
matter mass iMandelbaum et al.1 [iooebl '). The former are 
subject to systematic uncertainties due to the initial mass 
function (IMF) at the ~ 20-30 per cent level, whereas the 
latter is subject to ~ 10 per cent uncertainty in the mod- 
elling in addition to ~ 10-15 per cent statistical error. 

The selection of mass scales begins with the ansatz that 
we would like to approximately bracket the typical lumi- 
nosity and mass range of observed lensing systems. Given 
the range of lens galaxy luminosities in the CASTLES and 
SLAGS strong lens samples, we choose models with total 
luminosities of ~ 2I/« and ~ 7L* in the (SDSS) r-band. 
We refer to the lower-luminosity model interchangeably as 
"galaxy scale" or "lower mass scale" model. The higher- 
luminosity model is essentially a massive galaxy with con- 
vergence that is significantly boosted by a group dark matter 
halo, but we often refer to it as a "group scale" or "higher 
mass scale" model for brevity. In both cases, the model 
includes concentric stellar and dark matter components; 
thus, our current results cannot be applied to satellites in 
groups/clusters, which will be the subject of future work. 
We also neglect contributions of nearby (satellite) galax- 
ies to the surface density, focusing only on the host galaxy, 
and we neglect c ontributions from other structures along the 
line of sig ht (e.g.' Keeton fc Zabludo3l2004l : rMomcheva et al.1 
l2006l: IWi lliams c t al.1 120061). T he masses are selected from 
the lMandelbaum et al.l ( 2006bl ) SDSS weak lensing analysis 
using results for early-type galaxies as a function of lumi- 
nosity, as shown in Table 3 or Fig. 4 of that paper. 

The model parameters for these mass scales are sum- 
marised in Table [1] We define the virial radius entirely using 
comoving quantities, with riso satisfying the relation 



the universe (using Q,m ~ 0.27). The definition is in princi- 
ple arbitrary, but this choice was used to analyse the weak 
lensing results. In the weak lensing results, the profiles were 
found to be consistent wit h NFW for early-t ype galaxies on 
all scales studied (see also iMandelbaum et al, 2006a. which 
focuses on z ~ 0.25 Luminous Red Galaxies). Thus, not only 
our use of the masses, but also the form of the density pro- 
files, is observationally motivated outside of ~ 20h~^ kpc. 



4.2 Length scales 

Wo choose the dark matter halo concentration from the 



'Bullock ct al. (2001) result that Cdn 



10(M/M„i) 



-0.13 



af- 



180p = 



3Afi8o 



(31) 



where Migo is the virial mass and p is the mean density of 



ter converting to our mass definition and using inner slope 
7dm — 1, zl — 0.3, as = 0.75 and flm = 0.27. The nonlinear 
mass A/ni is the mass in a sphere within which the rms linear 
density fluctuation is equal to Sc, the overdensity threshold 
for spherical collapse. We obtain Cdm — 8.4 for the lower 
mass scale and Cdm — 5.6 for the higher mass scale (see also 
Table [1} . The scale radius of the dark matter component 

then follows as rs,dm = ri80,dm/cdm. 

We determine the scale length of the stellar components 
based on SDSS photometry of galaxies of the appropriate 
luminosity. Fitting de Vaucouleurs profiles, i.e., Sersic pro- 
files with index n = 4, to the growth curve of these elliptical 
galaxies yields the half-light (or effective) radius Re- Assum- 
ing a constant stellar mass-to-light ratio. Re is equal to the 
projected stellar half-mass radius R'f^ ^ as given in Table [1] 
For a Hernquist cusped density profile with (7*,n*) = (1,4) 
as we adopt for the stellar component, the scale radius fol- 
lows as rs,* — R'f^^/1.8153. 

Given the large scatter in dark matter halo concen tra- 
tion values in simulations (0.15 dex: lBullock et al]l200lh . we 
will also redo a subset of the analysis using concentrations 
lower and higher than the fiducial values used for the main 
analysis. This work will allow us to quantify biases associ- 
ated with halo concentration. Realistically, we expect the 
changes in the concentration Cdm to be somewhat degen- 
erate with changes in the inner slope 7dm, as both change 
the amount of dark matter in the inner parts, although there 
may be some differences because changing 7dm affects the in- 
ner logarithmic slope while changing Cdm does not. In Paper 
II we quantify the approximate degeneracy between changes 
in Cdm and 7dm for the dark matter component, while keep- 
ing the stellar component fixed. 

4.3 Density profiles 

For our main set of simulations, we use Hernquist and gener- 
alised NFW profiles for stellar and dark matter components, 
respectively (cf. Sections 12.31 and 12. 4p . While there is some 
uncertainty over the true dark matter profile in A''-body sim- 
ulations (particularly in the inner regions; see Section [2.5p . 
and the initial dark matter profile may in any case be mod- 
ified due to the presence of baryons, we use the generalised 
NFW model because it allows freedom in the form of the 
cusp. As discussed in more detail below (Section I4.5|l . we 
also examine how our results depend on the form of the 
density proflles by running a set of simulations that use de- 
projected Sersic proflles for the intrinsic (3d) densities of 
both the stellar and dark matter components. Here we fo- 
cus on the main set of simulations with cusped proflles. 
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Table 1. Summary of the spherical density profile parameters for galaxy and group scale models. With the exception of R'^ ^, all numbers 
given are intrinsic {3d). 



Description 


Symbol 


Galaxy-scale 


Group-scale 


Unit 


Dark matter component 


Virial mass 


-'^^180, dm 


3.4 X 10l2 


6.7 X 10^^ 






A-flSO.dm 


4.7 X 10i2 


9.3 X 10^^ 


Mq 


Virial radius 


''180, dm 


390.7 


1055 


comoving kpc 




''180, dm 


417 


1127 


physical kpc 




''180, dm 


97.1 


262 


arcsec 


Concentration 


Cdm (7dm = 1) 


8.4 


5.6 




Scale radius 


rB,dm (7dm = 1) 


49.7 


201 


physical kpc 




rs, dm (7dm = 1) 


11.56 


46.8 


arcsec 


Stellar component 


Total mass 


Aftot,* 


1.16 X 10" 


5.6 X 10" 


Mo 


Half-mass radius 




4.07 


11.7 


kpc 


(projected) 




0.95 


2.7 


arcsec 


Scale radius 


rs,* (7* = 1) 


2.24 


6.4 


physical kpc 




rs,* (7* = 1) 


0.52 


1.5 


arcsec 



4-3.1 Inner and outer slope 

For our main set of simulations with cusped density profiles, 
we fix tiie stellar component to have a Hernquist profile, with 
inner slope 7* = 1 and outer slope n* = 4. Based on cos- 
mological A^'-body simulations, we fix the outer slo pe of the 
dark matter simulation to ridm = 3 (e.g., Na varro &: et al.l 
I2OO4I ). 

The inner slope of the dark matter profile is more prob- 
lematic. While A'^-body simul ations can suKgest for ms for 
dark matter halo profiles (e.g.. lNavarro fc et al.ll2004 ). there 
is some question whether those profiles are affected by the 
presence of a stellar component, and if so, in what way. 
One formalism for analysing how a stellar component af- 
fects a dark matter profile is called adiaba t ic contraction 
(AC: lYound [l980l: iBlumenthal et al.l [l986l : iGnedin et all 
I2OO4I : [ Sellwood fc McGaughll2005l ). The physical picture is 
that as the baryons cool and condense into the centre of 
the system, they draw some of the dark matter in as well, 
causing the dark matter halo profile to steepen. Persuasive 
observational evidence for or against AC in the galaxies we 
are modelling does not yet exist, though there are claims 
that the theoretic al assumptions beh ind AC are not valid 
for these galaxies. iNaab et al.l (|2007l ) suggest that the op- 
posing effect, growth through dissipationless infall of matter 
onto the galaxy (which tends to push the dark matter out- 
wards from the centre) may cancel out the steepening due 
to baryonic cooling. We expect the effects of infalling mat- 
ter to be even more significant for our higher mass model, 
which is equivalent to a fairly massive group. Furthermore, 
hydrodynamic simulations that demonstrate that adiabatic 
contraction dominates may suffer from over-cooling of the 
baryonic component, so we cannot use them to determine 
the extent of the effect either. Thus, we cannot be certain 
how much the dark matter profiles may be modified due to 
the different processes involved in galaxy formation. 

We consider the range of possibilities in two ways. First, 
for the main set of simulations with cusped density pro- 
files, we examine three values of the dark matter inner slope: 



7dm = {0.5, 1.0, 1.5}. As we shall see, this allows for a broad 
range of lensing properties. Second, for a subset of simu- 
lations with 7dm ~ 1 NFW dark matter profiles, we ex- 
plicitly include the effec ts of AC, using the form alisms in 
IBlumenthal et al.l (|l986D and lCnedin et al.l (|2004l ). 



4-3.2 Normalisation 

When the inner slope of the dark matter 7dm is changed, the 
virial mass Miso also changes. To preserve M180 we must 
change either the density amplitude po or the scale radius 
rs of the dark matter (or both at the same time) . Preserv- 
ing in addition the (reduced) shear at the virial radius r-i8o, 
would break this "degeneracy" in the normalisation. Ifow- 
ever, for individual strong lens galaxies, weak lensing and 
other constraints on the density are very weak if present 
at all. Nevertheless, the different choices for the normalisa- 
tion might not result in significantly different lensing cross 
sections. 

To investigate this, we compute the change in total 
lensing cross-section with respect to the fiducial case with 
7dm ~ 1 when varying po and/or rs to preserve Misc. The 
results are summarised in Fig.[T]for both the lower (top) and 
higher (bottom) mass scale. The curves in the left panels 
show the relative difference in unbiased cross-sections, while 
the curves in the right panels are for biased cross-sections 
with a limiting luminosity of 4Lt (see Section 15.41 below) . 
The solid curves show the relative difference in cross-section 
when both rs and po axe kept fixed when changing the inner 
slope away from the fiducial value 7dm = 1, i-e., not pre- 
serving the virial mass A/iso. The dotted and dashed curve 
show the difference when preserving Miso by varying either 
Po (fixing rs) or (fixing po), respectively. 

The results for the unbiased and biased cross-section are 
similar. At 7dm ~ 0.5, the offsets between the solid, dotted 
and dashed curves corresponding to the three different nor- 
malisations are marginal, which is expected since for a shal- 
lower slope the contribution of the dark matter component 
decreases with respect to the stellar component which we 
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Galaxy scale, unbiased 
Fix r^, p„ 

Fix r^, M|ao - 

Fix p„, M,j,„ /<j 


Galaxy scale, biased / 

(L/U = 4, total)//* 

/ '' - 
/ '' ■' 

— //.-'_ 

/ /-■' 
/// 

A' 

r 


1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
I Group scale, unbiased i 

^^^^^^^^^^^^^ 

, 1 , , , 1 , , , 1 , , , 1 , , , 1 , 


1 1 1 1 1 1 1 1 1 1 1 1 1 1 

~ Group scale biased / 
; (L/L.='4, total)/^>; 

, , , 1 , , , 1 , , , 1 , , , 1 , r 



0.6 0.8 1 1.2 1.4 0.6 0.8 1 1.2 1.4 



Figure 1. Variation in total lensing cross-section for different 
normalisations for the lower, galaxy (top) and higher, group (bot- 
tom) mass scale. The change in cross-section is shown relative to 
that of the fiducial case with dark matter inner slope of 7;jni = 1- 
The panels on the left show the relative difference in unbiased 
cross-sections, while the panels on the right are for biased cross- 
sections with a limiting luminosity of 4L* (see Section 15.41 for 
details). The solid curves show the relative difference in cross- 
section when both Ts and po are kept fixed, while the dotted 
and dashed curve show the difference when preserving Afiso by 
varying po (fixing rs) or Vs (fixing po), respectively. 

keep fixed. However, the dark matter contribution becomes 
important when 7dm = 1-5, in particular for the high-mass 
scale as can be seen from the different scaling of the vertical 
axes. We see that in aU cases the solid curve leads to higher 
cross-sections than the dotted and dashed curve when we 
preserve Migo, which is expected given the increase in mass 
of tens of per cent. 

The latter two normalisations result in changes in the 
cross-section of which the relative difference is well within 
the combined statistical and systematic uncertainty in the 
cross-sections (i.e., while we have attempted to make these 
models fairly realistic, there are some systematic uncertain- 
ties, so we do not trust them to per cent level precision). 
This implies that varying po and/or Ts to preserve Mi8o 
does not change the lensing properties in a way that is 
significantly different. In what follows we choose to fix 
(and hence the concentration riso/^s), because it is sim- 
pler and in line with preserving the ( weak) concentration - 
mass relation for dark matter halos (e.g. iBuUock et al.ll200ll ) 
and the size-lum i nosity relation for ea rly-type galaxies (e.g. 
IShen et al.|[2003l : iBernardi et"al]|200'i1 ). as discussed further 
in Section [4.5.21 below. 

^.3.3 Intrinsic profiles 

For the spherical density profiles, we present a set of plots 
showing various quantities for both mass scales. We begin 



with Fig. [2l which shows the intrinsic, 3d mass density pro- 
file p{r) for the cusped density profiles (NFW plus Hern- 
quist), for both mass scales. It is apparent that for 7dm = 0.5 
and 7dm ~ 1 density profiles (stellar plus dark matter), the 
dark matter component is negligible for scales below ~ 2 
kpc for both mass scales, whereas the 7dm ~ 1-5 models 
have significant contributions from dark matter for all scales 
shown. The difference between these cases is somewhat re- 
duced when considering the profiles in projection along the 
fine-of-sight, due to the dominant dark matter contributions 
at large intrinsic radius. Nonetheless, we expect a significant 
break in the strong lensing properties going from 7dm = 1 
to 7dm = 1.5. 

The colour curves in Fig. [5] show the effects of AC on 
the cusped, 7dm ~ 1 density profiles. We can see (bottom- 
left panel) that for the lower mass scale, which has a higher 
ratio of stellar to halo mass, AC tends to increase the den- 
sity on scales below ~ 10 kpc by ~ 40 per cent depending 
on the model used, at the expense of a slightly decreased 
density (~ 3 per cent) for larger scales. For the higher mass 
scale (bottom-right panel), the de nsity below 10 kpc t ends 
to increase by ~ 30 per cent for the lCnedin et all (|2004l ) AC 
prescription, again at the expense of a percent-level decrease 
on much larger scales. 

Fig. |3] shows the negative logarithmic slope y{r) of the 
intrinsic mass density for both mass scales. We can see that 
for the individual components of the models, the logarith- 
mic slopes behave smoothly according to equation (|12[) . but 
for the full model, there can be complex, non-monotonic be- 
haviour. The latter "wiggle" in Fig. [3] is stronger for the 
cusped than the Sersic density profiles, refiecting the break 
in the cusped density profile while the Sersic density profile 
turns over smoothly. 

The intrinsic profiles of our galaxy models are consistent 
with the total mass density profile derived bv lCavazzi et al.] 
(2007), based on a joint strong-lensing and weak- lensing 
analysis of 22 early-type galaxies from the SLACS sample. 
The inferred intrinsic density profile (their Fig. 8b) has an 
amplitude which falls in between our lower-mass galaxy and 
higher-mass group scale (for an adopted Hubble constant 
h — 0.72). Moreover, the corresponding slope matches very 
well the wiggle in the logarithmic slope of our models in 
Fig. [3l Since this non-monotonic variation is around a value 
of 7 = 2 (indicated by the dashed horizontal line in Fig. O, 
our models as well as the total mass density inferred by 
ICavazzi et al.l (,2007 ) are rather close to isothermal over a 
large radial range around the Einstein radius. However, the 
density overall is not consistent with isothermal, with clear 
deviations at small and large radii, whi ch happen well within 
the two-decade radial range for which 'Ga vazzi et al.l (|2007| ) 
claim consistency with isothermal. We discuss this "isother- 
mal conspiracy" further when studying the profiles in pro- 
jection next. 

4.3.4 Projected profiles 

In addition to the plots of the intrinsic, 3d profiles, we 
also show projected, 2d quantities, which are observation- 
ally more relevant. Fig. |4] shows the surface mass density E 
as a function of the projected radius R' , in physical units 
(kpc) and at the top in angular units (arcsec) for a lens 
galaxy at the fiducial redshift of zl ~ 0.3. The total (dark 
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Figure 2. Intrinsic mass density versus intrinsic radius for cuspod (top) and deprojected Sorsic (middle) density profiles, for the lower 
(left) and higher (right) mass scale. Top panels: The dotted curve is the Hcrnquist profile of the stars, while the single-dot-dashed, 
solid and triple-dot-dashed curves represent the generalised NFW dark matter profiles with inner slopes 7(j,n of 0.5, 1.0 (NFW) and 1.5, 
respectively. The corresponding thick curves show the total mass density profile by combining that of the stars and dark matter. The blue 
and red curves show the total mass density profile (increased by a factor ten for clarity) taking into account adia batic contraction (AC) o f 
the fiducial cusped (NFW ) dark matter profile (shown as black curve for comparison) using the prescriptions of lBlumenthal et al.l lll986l ) 
and iGnedin et all l|2004l ) , respectively. Middle panels: The thin curves indicate the stellar and dark matter profiles which combined yield 
the thick curves. The solid curves are based on the fiducial Sersic indices for both components, whereas the two other curves indicate 
the variation in the mass density when adopting the limits in n^m a-nd ''vr (see Section I4.5|l . Bottom panels: The black curves show 
the relative difference in the cusped and deprojected Sersic mass densities for the stellar component (dotted), dark matter component 
(thin solid) and the combination (thick solid), using the fiducial density profile parameters. The blue and red curves show the relative 
difference between fiducial cusped dark matter profile and those after taking into account AC, with corresponding axis labels on the 
right side. The dotted vertical line indicates the scale radius Vs^i, of the stellar component, whereas the innermost and outermost solid 
vertical lines show respectively the scale radius S'lid virial radius rigg dm of the dark matter component, using the fiducial density 

profile parameters. 
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0.1 1.0 10.0 100.0 1000.0 0.1 1.0 10.0 100.0 1000.0 

r [kpc] r [kpc] 

Figure 3. The negative logarithmic slope of the intrinsic mass density as a function of the intrinsic radius, for cusped (top) and 
deprojected Sersic (middle) density profiles, and the relative difference (bottom), for the lower (left) and higher (right) mass scale. The 
horizontal dashed line indicates the slope of an isothermal density profile. The meaning of the other curves is the same as in Fig. [2] 



matter plus stars) for 7dm — 0.5 and 7dm = 1-0 are 

nearly identical and star-dominated for < 1 kpc. As antic- 
ipated, this is a smaller transverse scale than that for which 
they are identical in 3d. 

The colour curves in Fig.|4]show the effect of AC on the 
projected cusped, 7dm = 1 density profiles. As expected, the 
projection effects have brought in larger scales where AC is 
less significant, so the projected profiles with AC tend to be 
higher by ~ 30 (lower mass scale) and ~ 25 (higher mass 
scale) per cent on the scales of interest (a few kpc for the 
lower mass scale, ~ 10 kpc for the higher mass scale). 

Fig. [5] shows the negative logarithmic slope of the sur- 
face mass density, 7'(i?'). If the density profile were a pure 
power-law, the 2d and 3d negative logarithmic slopes would 
be related by 7' = 7 — 1. Comparing Fig. [5] to Fig. |3] we 



see that this simple relation does hold at small and large 
radii for our cusped models, since these models are asymp- 
totically power-laws. However, it is not valid at intermediate 
radii, because our models — especially the full models con- 
taining both stellar and dark matter components — are far 
from being simple power-laws. This is particularly notable 
in the vicinity of the Einstein radius, where the 3d slope 
varies much stronger with radius than the 2d slope. 

Fig. |6] shows the projected dark matter fraction /dm 
within projected radius R' . As shown, by varying the dark 
matter inner slope in the cusped models, we are able to 
explore models with widely varying fractions of dark mat- 
ter. The dark matter fractions in our models are consistent 
with the range of values inferred from attempts to com- 
pare stellar and dark matter masses in observed lens systems 
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Figure 4. Projected (or surface) mass density versus projected radius, for cusped (top) and deprojected Sorsic (middle) density profiles, 
and the relative difference (bottom), for the lower (left) and higher (right) mass scale. At the top, projected radius is given in angular 
units for a (lens) galaxy at the fiducial redshift of = 0.3, corresponding to an angular diameter distance of Dj^ ~ 900 Mpc. With a 
source at redshift zg = 2.0, the vertical dashed lines indicate the Einstein radii Rein for the composite cusped and deprojected Sersic 
models, using the fiducial density profile parameters. The dotted and solid vertical line indicate the half-mass radii ijj^ ^ and ij^ 
containing half of the total mass of the stellar component and half of the virial mass of the dark matter component, respectively. In case 
of a constant mass-to-light ratio, ijj^ ^^^^ is equal to the effective radius Re of the surface brightness. The meaning of the curves is the 
same as in Fig. [2] 
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Figure 5. The negative logarithmic slope of the surface mass density as a function of the projected radius, for cusped (top) and 
deprojected Sersic (middle) density profiles, and the relative difference (bottom), for the lower (left) and higher (right) mass scale. The 
horizontal dashed line indicates the slope of a (projected) isothermal density profile. The meaning of the other curves is the same as in 
Fig.H 



dRusin et all 1 2003': 'Ru sin fc Kochaneklliooa : iFerreras et aD 
l2005l : lBolton et al. 2008bl ). 

Finally, Fig. [7] shows the deflection angle a{R') = 
d(j>/dR' as function of the projected radius R' . The Einstein 
radii follow from the solution of a{Rcin) = ^cin, which cor- 
responds to the intersections of the dashed one-to-one line 
with the thick curves, indicated by the filled circles. The 
corresponding Einstein radii appear at about the half-mass 
radius R'f^ ^ , which in case of a constant stellar mass-to-light 
ratio corresponds to the effective radius Re- 

Figs. [5] and [7] contain an important point related to 
attempts to combine strong lensing and stellar dynamics 
to measure the density profile in the vicinity of the Ein- 
stein radius (e.g.. lTreu fc Koopmans|[2004l : iKoopmans et alj 



12006*). The slope of the composite surface mass density for 
models with 7dm of 0.5, 1.0 and 1.5 are all close to unity 
near the Einstein radius. For the lower-mass scale models 
(top-left panel) there is a little more spread around a some- 
what higher value 7' ~ 1.2 ± 0.2 than for the higher-mass 
scale (top-right panel) with 7' ~ 0.9 ± 0.1. However, all are 
close to the slope 7' = 1 of a projected isothermal den- 
sity as indicated by the horizontal dashed line in Fig. [5] In 
Fig. [T] the deflection angle for the higher-mass scale (top- 
right panel) increases slowly with radius around the Eintein 
radius, but for the lower-mass scale the curves are nearly 
flat around (and beyond) Rein, similar to a constant deflec- 
tion angle in case of an isothermal density. Even though 
the stellar density and dark matter density with different 
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Figure 6. The fraction of projected dark matter within a projected radius, for cusped (top) and deprojected Sersic (middle) density 
profiles, and the relative difference (bottom), for the lower (left) and higher (right) mass scale. The meaning of the curves is the same as 
in Fig. |4] 



inner slopes are clearly non-isothermal, the combined den- 
sity has a projected slope and corresponding deflection an- 
gle consistent with isothermal around the Einstein radius. 
This finding implies that it will be challenging to use strong 
lensing to constrain the dark matter inner slope 7dm for a 
two-component model without bringing in information from 
significantly different scales. 

We note that the apparent isothermality of our com- 
posite models is a completely unintended consequence of 
the mass models we have chosen, not something we have 
imposed by design. However, it provides support for the 
realism of our models, when compared against observa- 
tions. There is much evidence from strong lensing (especially 
in combination with kinematics) that early-type galaxies 
have total density profiles that are approximately isother- 



mal i n the vicinity of the Einste i n radius (e.g., Rusin et al.l 
20031: iTreu fc Koopman3 l2004l : iRusin fc Kochanekl booa 



Koopmans et al.l l2006l ). The fact that a superposition of 



non-isothermal stellar and dark matter profiles can pro- 
duce a total profile that is approximately isothermal over 
a wide range of scales — known as the "isothermal con- 
spiracy" — is related to the tendency towards fiat rotation 
curves originally used to argue for a dark matter compo- 
nent. Finally, the same set of SLACS lenses can be fit with 
a sin gle power-law profile that is consistent with isother- 
mal (|Koopmans et alll2006l ). or by a more physically intu- 
itive combin ation of a stellar Hernqu ist profile and a dark 
matter halo l| Jiang fc KochanekllioOTh . All in all, the quasi- 
isothermal nature of our models appears to be consistent 
with interpretations of strong lensing by real galaxies. 
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Figure 7. Deflection angle (or derivative of the projected potential) as function of projected radius, for cusped (top) and deprojected 
Sersic (middle) density profiles, and the relative difference (bottom), for the lower (left) and higher (right) mass scale. The intersections 
of the diagonal one-to-one (dashed) line with the (thick) curves of the combined stellar and dark matter components yield the Einstein 
radii Rein as indicated by the filled circles and the corresponding vertical lines. The meaning of the other curves is the same as in Fig. |4] 
The deflection curves of the stellar and dark matter component separately show respectively a strong decrease and increase towards 
larger radii. However, the combination results for the lower mass scale in an almost flat deflection curve around the Einstein radius and 
beyond, similar to a constant deflection angle in case of an isothermal density profile. 



4.4 Density shape models 



We consider seven different models for the shape of tiie 
galaxy density: spherical, moderately oblate, very oblate, 
moderately prolate, very prolate, triaxial, and a mixed 
model described below. In all but the mixed model, the 
axis ratios used for the dark matter and stellar components 
are related in a simple way (see the full description in Sec- 
tion |4]4]T]). In all cases the axes of the two components are 
intrinsically aligned in three dimensions. 



4-4- 1 Axis ratios 

We adopt the mean values (6/a)dm = 0.71 and (c/a)dm = 
0.50 from jing fc Suto (,2002 ) for the axis ratios of the triax- 
ial density of the dark matter component. These ratios were 
determined from cosmological A'^-body simulations, without 
any stellar component. We set (6/a)dm = 1 for the oblate 
shapes, with the same (c/a)dm = 0.50 for the moderately 
oblate case, and (c/a)dm = 0.71 x 0.50 ~ 0.36 for the 
very oblate case. We set (fe/a)dm = (c/a)dm for the pro- 
late shapes, with (c/a)dm = 0.71 for the moderately prolate 
case and (c/a)dm = V0.71 x 0.51 ~ 0.60 for the very pro- 
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late case. Note that in the very oblate and very prolate case 
these choices preserve the product (6/a)dm x (c/a)dm of the 
triaxial case, and hence the enclosed ellipsoidal mass. 

Next, we use results from the cosmological gasdynamics 
simulations in iKazantzidis et al] (|2004l ) to make a (crude) 
conversion from the dark matter shape to the rounder shape 
of the stars: {b/a)t = 0.6 (fe/a)dm + 0.4, and similarly for 
{c/a)i,. In practice, in that paper, the axis ratios were a 
function of separation from the centre of the halo; we neglect 
those effects in this work. Also, in that work the dark matter 
halos were themselves made more round due to the presence 
of baryons, with up to ~ 15 per cent effects out to the virial 
radius. However, those simulations may suffer from baryonic 
overcooling, which would tend to accentuate these effects, 
so we do not incorporate the rounding of the dark matter 
component at this time. 

The final model, the mixed shape model, is constructed 
using the very oblate stellar shape combined with the triax- 
ial dark matter halo shape. Part of t he motivation f or do - 
ing so is based on the results given bv lLambas et all (Il992l ) 
for the axis ratio distributions and 3d shapes of early-type 
galaxies as inferred from the distributions of their projected 
shapes. In that paper, the elliptical sample containing 2 135 
galaxies was found to have a projected shape distribution 
consistent with a (weakly) triaxial intrinsic shape distribu- 
tion. Fitting Gaussian distributions in both axis ratios, they 
found a best-fitting mean and standard deviation of 0.95 and 
0.35 in (&/a)*, and 0.55 and 0.20 in (c/a)*. The 4 782 SO 
galaxies in their sample had a projected shape distribution 
with best-fit mean value (fe/a)* = 1 consistent with oblate 
axisymmetry, and with best-fitting mean and standard de- 
viation of 0.59 and 0.24 in (c/a)*. 

Further motivation is provided by detailed dynami- 
cal models of nearby elliptical and lenticular galaxies that 
accurately fit their observed surface brightness and (two- 
dimensional) stellar kinematics. These models show that 
the inner parts of lenticular as well as many elliptical 
galax ies, collectively called fast-rotators (lEmsellem fc et al] 
l2007l). are consistent wit h an oblate axisymmetric shape 
i Cappellari fc et al.|[2007l ). The slow-rotator ellipticals, on 
the other hand, might be rather close to oblate axisymmet- 
ric in the centre but rapidly become t ruly triaxial at larger 
radu (|van den Bosch et al.ll2008l . |2008| ). 

We thus consider this mixed model, with a triaxial dark 
matter shape as in cosmological A'^-body simulations plus an 
oblate stellar shape as inferred from observations, as being 
potentially consistent with reality for the lower mass scale. 
The triaxial shape for both dark matter and stars may be 
more appropriate for the higher mass scale. 

The resulting axis ratios for each shape model, both for 
the dark matter and stellar component, are given in Tabled 
While we do not explore a broad range of axis ratios due to 
the prohibitive numbers of simulations involved, we can still 
get some sense of the trend with axis ratio by looking at 
the changes from spherical to moderately and very oblate or 
prolate shapes. 



Scale lengths for non-spherical density profiles 

When constructing the non-spherical density profiles, we fix 
the concentration to be the same as in the spherical case, 
but set the scale length a such that the mass within the 



Table 2. Adopted axis ratios of the ellipsoidal intrinsic mass den- 
sity for the dark matter and stellar component for seven different 
shape models. 



Shape 


(Va)dm 


(c/a)dm 






Spherical 


1.00 


1.00 


1.00 


1.00 


Moderately oblate 


1.00 


0.50 


1.00 


0.70 


Very oblate 


1.00 


0.36 


1.00 


0.62 


Moderately prolate 


0.71 


0.71 


0.83 


0.83 


Very prolate 


0.60 


0.60 


0.76 


0.76 


Triaxial 


0.71 


0.50 


0.83 


0.70 


Mixed 


0.71 


0.50 


1.00 


0.62 



ellipsoidal virial radius miso = riso/a, is equal to the virial 
mass Afi8o. This is done using 



{b/a) (c/a) = r'i 



(32) 



In this case, the length scale of the halo is rescaled while 
maintaining its concentration riso/rs- 



4.4-3 (Mis) alignment of projected axes 

We assume that the intrinsic axes of the stellar and dark 
matter component are aligned. For the axisymmetric models 
this means that the minor and major axes of the surface 
mass densities of the stars and dark matter are aligned as 
well. However, in case of a triaxial dark matter component 
the projected axes are misaligned with respect to those of 
the intrinsically rounder stellar component. 

The misalignment Atp = IV'dm — follows from equa- 
tion © for the different intrinsic axis ratios of the dark 
matter and stellar component (Table [2} as a function of the 
polar and azimuthal viewing angles (??, ip). The left panels of 
Fig.[5]show the resulting Atp for the triaxial (top) and mixed 
(bottom) shape model. Depending on the viewing direction, 
the projected axis ratios of the stars and dark matter can be- 
come completely orthogonal {Atp = 90°). At the same time, 
however, the projected ellipticity of the stars 1— (&'/a')t goes 
to zero, as shown in the right panels of Fig. [S] This means 
that the stars appear round on the plane of the sky, so that 
their position angle on the sky is very hard to establish. 

Even ignoring the stellar ellipticity, the probability of a 
significant misalignment is rather small in the case of ran- 
dom viewing direction^ For the triaxial shape model only 
~ 1.2% of the area on the viewing sphere leads to misalign- 
ments Alp ^ 10 ° . This number increases to ~ 50.8% for the 
mixed shape model, but at the same time the area on the 
viewing sphere that leads to a round appearance of the stars 
is also larger: ~ 8.3% of the area has (&'/a')* ^ 0.95, while 
this is only ~ 2.1% for the triaxial shape model. 

Our use of models that lead to almost no situations 
with significant projected ellipticity and misalignment of 
DM and stellar components is consistent with studies of 
relatively isolated lens systems, for which the position an- 
gles of the light a nd of the projected total mas s agree 
to within < 10° (|Keeton et al ] Il998l : iKochanekl l2002al : 



^ We show in Paper II that there is a selection bias in orientation 
for non-spherical lens galaxies. 
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Figure 8. The left panels show the misalignment Aip between the projected axes of the surface mass density of the stellar and dark 
matter component of the triaxial (top) and mixed (bottom) shape models, as a function of the polar viewing angle i9 and azimuthal 
viewing angle (p. In practice, Aip is estimated from modelling of the lensing images and the position angle of the stars on the sky. The 
latter becomes very uncertain if the stars appear too round on the plane of the sky, i.e., in case of a very small ellipticity 1 — (6'/a')* 
shown in the right panels. Note that these are not equal-area plots since the horizontal axis is rather than cos 6. 



iKoopmans et al.ll2006l : iBolton et al]|2008lj ) . Moreover, both 
the observed surface brightness and stellar kinematics of 
early-type galaxies are well-fitted by dynamical models with 
the intrinsic density parametrised by multiple (Gaussian) 
com ponents that have diff e rent shapes but are all c o-axial 
(e.g. lEmsellem eral]|l994l : ICappellari fc et aLllioO^ . This 
includes the photometric and ki nematic twists observed in 
parti cularly the giant ellipticals (|van den Bosch et al.|[200^ . 
l2008h . 

The above is all empirical evidence in favour of an ini- 
tial guess involving a high degree of alignment. While there 
are many reasons why these assumptions are not perfect, 
here we do not explore the effects of intrinsic misalignment 
between the shapes, because of the enormous amount of pa- 
rameter space involved. Future work should explore this is- 
sue further, as well as the issue of warps in the intrinsic 
and/or projected profiles. 

4.5 Deprojected Sersic simulations 

In addition to the simulations with cusped density profiles, 
we also create a set of simulations with deprojected Sersic 
density profiles (Section 12. 5|) . which also provide a good de- 



scription of both stars and dark matter over a range of mass 
scales (Section I2.3|l . As before, we choose the different pa- 
rameters for the dark matter and stellar density to create 
observationally realistic galaxy models. 



4-5.1 Density parameters 

For the total mass of the stellar component we use the same 
values given in Table [l]for the lower and higher mass scales. 
Even though the total mass of the dark matter component 
is also finite in this case, we adopt for consistency the same 
virial mass scale by setting the mass within the virial radius 
'"180, dm equal to Migo,dm as given in Table [1] In spherical 
geometry, the scale radius of the Sersic density is the half- 
mass radius, which for the stellar component is ^ given in 
Table [T] This, in turn, is equal to the effective radius Re in 
the case of a constant stellar mass-to-light ratio. For the dark 
matter component we numerically derive the radius R'f^ 
that encloses a projected mass that is half of the virial mass 
Afi80,dm of the spherical NFW dark matter halo, and use the 
same half-mass radius for the deprojected Sersic profiles. 
This procedure results in R'h dm = 2.2699 r^, dm — 112 kpc 
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and -RJ, dm = 1.6882 Vs.dm — 339 kpc for the lower and higher 
mass scale, respectively 

The only parameter left to choose is the Sersic index 
n. To derive the Sersic indices n* of the stellar component, 
we use the relation logn* = -0.10 (M b + 18) + 0.39 for 
Virgo early- type galaxies as derived by iFerrarese fc et al] 
l|2006l ). A luminosity of L* corresponds to an absolute John- 
son S-band and SDSS r-band m agnitude of about M b* = 
-19.4 + Slog (/!. = 0.72) = -20.1 (iFaber fc et al.ll2007l) and 
Mr, = -20.28 -I- 5 log {h = 0.72) = -21.0 (|Blanton fc etlo] 
I2OO3I ). respectively. With B — r = 1.32 as a typical colour 
for elliptical galaxies l|Fukugita et al.lll995l ), this means that 
our choice for the lower and higher mass scales of ~ 2L, and 
~ 71/, in the r-band translate to ~ 1.31/* and ~ 4.71/, in 
the B-band, or about Mb = -20.4 and Mb = -21.8. We 
thus obtain Sersic indices n* — 4.29 and = 5.87 for the 
stellar component of the lower and higher mass scale, respec- 
tively. For the dark matter component, we use the Sersic fits 
to sim ulated ACDM halo profiles presented in lMerritt et al.l 
l|2005l) . Taking into account the difference in virial mass def- 
inition, we obtain nam = 2.96 and nam = 2.65 for the lower 
and higher mass scales, respectively. 

In addition to these fiducial Sersic indices, we also con- 
sider values consistent with the observed scatter. For a given 
a bsolute luminosity Mb, t he scatter in the observed logn^, 
in IFerrarese fc et al.l l|2006l ) is approximately ±0.1 dex. This 
scatter results in a range in n* of [3.41, 5.40] and [4.66, 7.39] 
for the lower and higher mass scales, respectively. Similarly, 
for a given vi r ial m ass Migo, the range of fitted log nam in 
iMerritt et al.l (|2005l ) is approximately covered by ±0.05 dex. 
Since these findings are based on a low number (19) of sim- 
ulations, we consider indices ±0.10 dex around the fiducial 
values. This scatter implies a range in nam of [2.35, 3.73] and 
[2.10, 3.33] for the lower and higher mass scales, respectively. 

4-5.2 Normalisation 



For the cusped density profiles in Section r4.3.2l changing the 
dark matter inner slope from the fiducial value 7dm = 1 only 
affects the inner parts of the density profile. As a result, the 
corresponding change in virial mass Migo,am (about —20% 
and +30% when changing to respectively 7am = 0.5 and 
7dm ~ 1.5), could easily be compensated by changing the 
density amplitude po and/or the scale radius Vs- For the de- 
projected Sersic density profiles, changing the index n affects 
the whole density profile rather than just the inner parts, so 
that the corresponding changes in the mass are much larger. 
Changing the fiducial dark matter index nam to the limits 
of the adopted range (while keeping fixed), results in a 
factor ~ 3 change in the virial mass Migo,am for both mass 
scales. Changing the fiducial stellar index to the limits of 
the adopted range, results in a factor > 100 change in the 
total stellar mass Mtot,<r. Preserving the mass thus requires 
very large variations in the central density po and/or the 
half-mass radius 7?Jj. 

Certainly for the stellar component, changing R'f^ is not 
an option, sin ce early-type galaxies obey a size-lumin osity 
relation (e.g. IShen et all l2003l : iBernardi et"aLl |20o3). Al- 
though there is some freedom due to scatter in this rela- 
tion and the stellar mass-to-light conversion, this relation 
excludes larger variations in -Rj, ^ when preserv ing Mtot.*. 
A similar conclusion is reached based on the iKormendvl 



(|l977l ) relation (h) cx R7°'^^, where (h) = Lt^t/iirRl) 
is the mean surface brightness within the effective radius 
Re- In case of spherical symmetry and a constant mass- 
to-light ratio R'f^ = Re and Ltot oc poRef{n*), where 
the latter dependence on the Sersic index follows from 
equation (|15p . Therefore, obeying the Kormendy relation 
means po ^i^^"^"'' ,/(«*) = cnst, while preserving the to- 
tal stellar mass implies po^e = cnst. Since q 



3 even at high(er) redshi f t (e.g . iLa Barbera et al] l2003l : 
Idi Serego Alighieri fc et ahl 120051 ), these two relations to- 
gether imply that R^ remains to be approximately constant, 
while Po can be varied to compensate for the change in n*. 

A similar scaling relation seems to hold for NFW dark 
matter halos between their concentration cam (or scale ra- 
dius in case of a constant virial radi us), and virial mas s 
A/i80,dm of the form cam oc M-°g^^^ jSullock et al.ll200ll ). 
Since the dependence on Migo.dm is weak and the scatter 
in the relation is rather large, we could not use it before 
in Section 14.3.21 when changing the inner slope 7am of the 
cusped density profile to break the degeneracy between the 
different normalisation options (which we showed in case 
of preserving A/igo,am was not necessary, since the effect 
of the diffe rent normalisations on the lensing cross-sections 
is similar). ISalvador-Sole et al.l (|2007l ) show that when an 
Einasto profile is used instead of a NFW profile, a similar 
concentration-mass relation exists for dark matter halos. Be- 
cause deprojected Sersic profiles fit simulated dark matter 
halos equally well as Einasto profiles (Section l2.3|l , we expect 
also in that case a correlation between cam and Migo.dm. As 
mentioned above, for the Sersic density the effect of chang- 
ing nam is much stronger than changing 7am for the cusped 
density, so that varying cam (while fixing po) to preserve 
A/"i8o,dm would violate this concentration-mass relation. As 
a result, for both the stellar and dark matter component we 
vary po to preserve the (total and virial) mass, while we keep 
the half-mass radius R'f^ fixed. 



4.5.3 Density profiles 

The profile of the spherical deprojected Sersic mass density 
p(r) is shown in Fig. [2] in the middle row and compared to 
the cusped mass density in the bottom row, for the lower 
and higher mass scales. The thin curves indicate the stel- 
lar and dark matter profiles which combined yield the thick 
curves. The solid curves are based on the fiducial Sersic in- 
dices for both components, whereas the two other curves 
indicate the variation in the mass density when adopting 
the above limits in Wdm and n-^. In a similar way, we show 
in Fig. Othe corresponding negative logarithmic slope ^{r), 
in Fig.|4]the Sersic surface mass density E(7?'), in Fig.[5]the 
corresponding negative logarithmic projected slope ')'(R'), 
in Fig.[S]the projected dark matter fraction /am(< ^-iid 
finally in Fig. [7] the deffection angle. 

As for the cusped density in Section 14.3.41 we find that 
both the slope of the composite surface mass density and the 
defiection angle (Fig. [S] and Fig. [7| are consistent with an 
isothermal density around the Einstein radius. Especially 
for the lower-mass scale, the projected slope for all cases 
converges to 7' = 1 around -Rein, and the curves of the de- 
flection angle are extremely flat at and beyond Rein. As for 
cusped models, this finding implies that constraints from 
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strong lensing cannot trivially be extrapolated to smaller 
and larger radii. 

4.6 Computation of surface mass density maps 

The lensing properties of a galaxy follow from its surface 
mass density T.{x',y'). For the above galaxy models this is 
the sum of the surface mass density of the stellar and dark 
matter components, each of which can be computed as fol- 
lows. For a given intrinsic shape (a, b, c) seen from a viewing 
direction 95), the projected semi-axis lengths {a',b') are 
given by equations (O-Q. At each position {x',y') on the 
plane of the sky, equation Q then yields the elliptic radius 
m' . The surface mass density follows by evaluation of equa- 
tion (O for a cusped density p(m) given in equation @, 
or in the case of a Sersic density directly as E(m') given in 
equation (|16|) . 

However, rather than doing the numerical evaluation 
for each position and viewing direction, we instead compute 
E(m') on a fine one-dimensional grid in m' , once for each in- 
trinsic density profile and shape. We then derive E(a:', y') by 
linearly interpolating onto this grid at the m' value that (ac- 
cording to equations|4HZl) corresponds to the position (x' ,y') 
and chosen viewing direction (tD, tp). In particular, in the case 
of the cusped density, this procedure means that we can 
avoid doing the numerical evaluation of equation (|3]) on a 
two-dimensional grid repeatedly for each viewing direction. 

Along with the construction of the cusped and Sersic 
galaxy models described above, we implemented and tested 
this efficient computation of the surface mass density in 
We use a linear grid in m' , except in the inner parts where we 
sample logarithmically in m! to accurately resolve a central 
peak in the density. Even so, when we compute the surface 
mass density on a regular square grid in {x' ,y') with a given 
dimension (box size) and sampling (resolution), the result- 
ing fixed pixel size might not be small enough to account 
for a strong central cusp. This is not the result of the in- 
terpolation but is due to the discreteness of the map-based 
approach. The effects of finite resolution and box size on 
the surface mass density and corresponding lensing proper- 
ties are further tested and discussed in Sections 15.71 and 15.81 

We choose the grid in {x',y') such that the I'-axis is 
aligned with the major axis of the surface mass density of 
the dark matter component. Even though we assume that 
the intrinsic axes of the stellar and dark matter compo- 
nent are aligned, in case of a triaxial dark matter com- 
ponent the projected axes are misaligned with respect to 
those of the intrinsically rounder stellar component. As dis- 
cussed in Section 14.4.31 above, the resulting misalignment 
angle Aip = |i/'dm — tp*\ follows from equation ((2]), and de- 
pends on the viewing direction {■&, tp). To take this misalign- 
ment of the stellar component with respect to the dark mat- 
ter component into account, we first rotate (x' ,y') over the 
misalignment angle At/; 

Xp = a;' cos At/) -I- y' sin A?/>, 
yp = — a;' sin At/) -I- cos Ai/i, 

and instead of equation Q, we then use m'^ = (xp/a')^ -|- 



http: / /www. ittvis.com/ProductServices/IDL.aspx 



(t/p/6')^ to define the elliptic radius m' . In case of the ax- 
isymmetric models Ai/) = 0, so that the stellar and dark 
matter component are always aligned. Finally, we add the 
maps of the stellar and dark matter components together 
to arrive at the total surface mass density T.{x' ,y') of the 
galaxy model. 

In Fig.[9l we show examples of derived surface mass den- 
sity maps for lower, galaxy-scale mass models with cusped 
density profiles. The figure includes, from top to bottom, 
the very oblate mass model viewed close to face-on and close 
to edge-on, and the mixed shape model. The latter is seen 
close to the viewing direction that results in the maximum 
combination of ellipticity of the stars (left column) and mis- 
alignment with respect to the dark matter (middle column) , 
resulting in a twist in the combined surface mass density 
(right column). In Section [5]9] we address the technical issue 
of how many viewing angles we need to sample in order to 
accurately characterise how the lensing behaviour changes 
with viewing direction. 



5 LENSING CALCULATIONS 

In this section, we describe in detail how we do the lensing 
calculations necessary for our analysis. We present numerical 
methods designed to balance efficiency with accuracy (Sec- 
tions [5Tl}{53}, together with extensive tests to validate those 
methods fSections l5.6H5.9]) . All of the calculations described 
here can be done with GRAVLENS version 1.99k, w hich is a 
version of the original software bv iKeetonI l|2001bl ) that has 
been extended for this projedQ. 

5.1 Lens potential, deflection, and magniflcation 

Given the two-dimensional convergence distributions for our 
models, we need to solve the Poisson equation (|22p to find 
the lens potential. In the present work, we actually need 
the first and second derivatives of the lens potential (for the 
deflection and magnification, respectively), and do not work 
with the potential directly. We use different approaches to 
obtain the derivatives, depending on the model properties. 
Note that we apply the methods discussed here to the stellar 
and dark matter components separately, and construct the 
composite lens model from a simple linear superposition of 
the mass components. 

If the model has circular symmetry and the conver- 
gence is known analytically, the deffection can be obtained 
from the one- dimensional integral over k{6) given in equa- 
tion (|25|) . Often this integral can be evaluated analytically, 
so the lens model is fully analytic. This procedure is used for 
circular deprojected Sersic models. (Fo r more detail a bout 
lensing by circu lar Sersic models, see ICardond [iooil and 
iBaltz et al.ll2009l .') 

If the model has elliptical symmetry and the conver- 
gence is known analytically, the first and second derivatives 

^ The features that are new in GRAVLENS vl.99k are identi- 
fied in the text. All the new features have been tested, and 
have instructions available through the GRAVLENS help com- 
mand, but are not fully documented in the manual. Version 
1.99k can be accessed using the link for latest updates at 
http : //redf ive . rutgers . edu/'keeton/gravlens. 
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Figure 9. Examples of surface mass density maps S as a function of position (x', y') on the sky-plane in case of lower, galaxy-scale mass 
models with cusped density profiles. From top to bottom, the very oblate mass model viewed close to face-on with inclination i ~ 14° 
and close to edge-on with i ~ 81°, and the mixed shape model viewed at a polar angle i? ~ 50° and azimuthal angle (/3 ~ 13°. The latter 
is close to the viewing direction that results in the maximum combination of ellipticity of the stars (left column) and misalignment with 
respect to the dark matter (middle column), resulting in a twist in the combined surface mass density (right column). The contour levels 
correspond to a factor of ten decrease in S going outwards, clearly showing that the stars are more concentrated than the dark matter. 



of the lens potential can be obtained from a set of one- 
dimensional integrals ove r ti(m' ) and its derivative dn/dm' 
(|Schramm|[l990 : Keetoiil liooTal ) . For elliptical deprojected 
Sersic models, we evaluate the integrals numer icallyI3yield- 



We note that the steep central profile of Sersic models can 
lead to numerical precision requirements that are more stringent 
than usual. With inadequate precision, there is a systematic and 
cumulative error that effectively changes the density profile. We 
achieve the necessary precision by decreasing the GRAVLENS inte- 
gration tolerance parameter inttol from its default value of 10~® 
down to 10"^''. 



ing what we call numerical lens models. Circular and ellip- 
tical Sersic models are new features in GRAVLENS that have 
been added for this project. 

For some models, not even the convergence is known 
analytically, so the projection integral in equation ([3)l must 
be evaluated numerically. In such cases we compute the con- 
vergence k{R') on a map and then solve the Poisson equa- 
tion using Fourier methods. Here, K = {x',y') is a two- 
dimensional angular position on the sky in the lens plane, 
with k' its counterpart in Fourier space. We use FFTs to 
find the Fourier transform K{k') of the convergence. We 
then rewrite the real-space Poisson equation ([22} in Fourier 
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space as 

^k'^F{k') = 2K{k') , (33) 

which is readily solved for the Fourier transform F{k') of 
the lens potential (j>{R'). Derivatives of the lens potential are 
easily computed in Fourier space: for example, the Fourier 
transform of V(j>{R') is ik'F{k'). We can then use inverse 
FFTs to convert back to real space. All of this analysis can 
be done with the GRAVLENS routine kap21ens (which has 
been added and tested for this project). We refer to lens 
models computed this way as map-based lens models, and 
we use them for our cusped models 

There are four important technical details to consider 
in the Fourier analysis. The first issue is avoiding "edge ef- 
fects" that may arise because we are using a rectangular box 
with an elliptical density map, and because FFTs naturally 
impose periodic boundary conditions. We embed the K-map 
in a grid at least twice as large as the desired box in each 
direction (and then round up to the next power of two, as 
required for FFTs). We take advantage of the fact that our 
density models have circular or elliptical symmetry, and we 
know (from Newton's third law) that the gravitational force 
inside an elliptical shell due to that shell is zero. We can 
therefore set the density to zero outside of a reference el- 
lipse that is larger than the region where we wish to study 
lensing, without affecting the lensing deflections or magnifi- 
cations (or differential time delays, although we do not study 
those). Specifically, if the box side length is L, and the angle 
of the projected mass distribution relative to the x' axis is tp 
(see also Section [2.21 and equation [2] in particular), we first 
define the axes of the projected mass distribution 

Xp = x' cos tp + y sin i/) 

yp = —X sin 'i> -\- y cos i/) 

and then set k to zero when 

<'''= + fe)' >(§)'■ '«> 

The second issue arises because equation (|33|l is ill- 
behaved at k' = 0. We avoid this problem by setting 
all Fourier transforms to zero at k! = 0. Since if(0) oc 
J k{R') dR' , putting K(0) = amounts to requiring that 
the total mass in the FFT box vanish. That, in turn, means 
the code effectively adds a negative, uniform mass sheet to 
offset the positive mass in our model K-map. We can com- 
pensate by explicitly inserting a positive mass sheet to com- 
pensate for the effective negative mass sheet. We have care- 
fully tested this approach to verify that it yields accurate 
lensing calculations (see Section [5. 6p . 

As a third issue, we need to make sure our results are not 
affected by the resolution of the grid on which we compute 
K. Finally, we must necessarily work in a finite box, and we 
need to check that our results are not affected by the size of 
the box. We examine the latter two issues carefully below in 
Sections 15.71 and 15.81 respectively. 

Hernquist and NFW models could be treated analytically in 
the spherical case and numerically in the elliptical case, since the 
convergence is known analytically. However, we must treat gen- 
eralised NFW models with 7 7^ 1 with the map-based approach, 
since the convergence is not analytic. We opt to handle all the 
cusped models in the same way for uniformity. 



5.2 Image configurations 

The GRAVLENS software features a general "tiling" algorithm 
to solve th e lens equation and find the images of a given 
source (see lKeetonll2001bl for details). We can check the re- 
covered images using general considerations from lens the- 
ory. First, we classify the images as minima, maxima, and 
saddle points as discussed in connection with equation (|27|l 
in Section [3. II We can then categorise the image configura- 
tion as follows: 

• A 2-image lens has one minimum and one saddle point 
(plus one maximum that we do not analyse). 

• A 3-image lens has two minima and one saddle point. 

• A 4-image lens has two minima and two saddle points 
(plus one maximum that we do not analyse). 

Note that all three cases satisfy the image number relation 
in equation (|28|l . Recall from Section 13.11 that we find but 
do not analyse maximum images because they are rarely 
observed and do not play a major role in most strong lensing 
applications. 

According to lens theory, these are the only viable 
multiply-imaged configurations that can be produced by an 
isolated galaxy whose central surf a ce density is shallowe r 
than E oc R'~^ (see iBurkd Il98ll : ISchneider et al.l \l99i ). 
Therefore, any set of recovered images that does not fall 
into one of these three categories must represent a numeri- 
cal error. 

In general, the error rate in GRAVLENS is very low. The 
few errors that do occur arise from two causes, which are 
related to the tiling algorithm. First, when the source lies 
very close to the caustic, two of the images can lie so close 
together (straddling the critical curve) that they are difficult 
for the tiling algorithm to resolve. Second, the tiling algo- 
rithm has adaptive resolution, which is valuable but presents 
certain challe nges in regions where the resolution changes 
(cf. Fig. 4 of i KeetonI [2001bh . In most cases checking the 
image configuration identifies the rare errors. There is one 
error that can slip past this check: if the code fails to find 
one saddle point of a 4-image lens, it will misclassify the 
system as a 3-image lens. We shall see below (Section I5.9|l 
that this error does occur, but at a low rate that does not 
corrupt our results at a significant level. 

5.3 Point source cross-sections 

Our initial focus will be on understanding intrinsic selec- 
tion biases in strong lensing surveys (see Paper II). We 
therefore study lensing cross-sections with the idea that 
any physical effect that enhances the cross-section will be 
more likely to appear in lens samples, while any effect that 
reduces the cross-section will be disfavoured. Historically, 
most lenses found in systematic surveys had point-like im- 
ages (quasars or radio sources). That is changing with th e 
advent of surveys hke SLAGS (|Bohon et all |2006| . l2008al ). 
but the statistics of galaxy/galaxy strong lensing are much 
mor e complicated than the statistics of point source lens- 
ing toobler et al]|2008l ). so we currently focus on the point 
source case. 

At the most basic level, the lensing cross-section is just 
the area of the source plane that leads to multiple imaging. 
We call this the total cross-section (crtot), and subdivide it 
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into the cross-sections for the regions that produce 2, 3, or 
4 images (cr2, (73, and 04). In real surveys, we also need to 
account for the fact that brighter lensed images are easier 
to observe; this "magnification bias" effect is discussed in 
Section 15.41 The cross-sections have dimensions of area on 
the sky, so the natural unit is arcsec^. 

We compute the cross-sections using Monte Carlo tech- 
niques. We place a source at a random position behind the 
lens, and solve the lens equation numerically to find all the 
images of this mock lens. We repeat this process many times, 
tabulate the number of mock 2-image, 3-image, and 4-image 
lenses, and use the input number density of sources to com- 
pute the associated cross-sections. We find that using 5000 
sources provides a good balance between run time and sta- 
tistical uncertainties. 

The most obvious approach is to distribute sources ran- 
domly and uniformly in the smallest circle enclosing the 
lensing caustics in the source plane. However, this approach 
yields poor sampling of high-magnification systems that are 
rare but important, especially for ma gnification bias (see 
Section[5]4|. For an alternate approach. iKeeton fc Zabludofj 
l|2004h pointed out that picking positions randomly and uni- 
formly in the image plane, and then mapping them to the 
source plane, yields a set of source positions that are ran- 
dom but weighted by the total magnification of all images 
(see Fig. nop . We must account for the magnification weight- 
ing in order to compute the cross-section and magnifica- 
tion bias properly, but that is a straightforward task. Com- 
pared with uniform source sampling, magnification-weighted 
source sampling is no more complicated in concept or algo- 
rithm, and yields superior results for Monte Carlo calcula- 
tions of lensing cross-sections and magnification bias. 

We determine the statistical errors on the cross-sections 
using bootstrap resampling. We create 500 resampled sets 
of mock lenses, and use the variance among them as the 
statistical uncertainty. All of this analysis (including boot- 
strap resampling) can be done with the GRAVLENS routine 
mockcsec (new in version 1.99k). 



5.4 Magnification bias 

Most lensing surveys have a flux limit that influences the 
statistical distribution of lenses. The fact that brighter im- 
ages are ea sier to observe lead s to a lensing "magniflcation 
bias" (e.g.. [Turner et al.lll984h . We account for this fact by 
computing not just the lensing cross-section itself, but the 
combination of cross-section and magnification bias, which 
we call the biased cross-section. 

Consider a survey in which the number of sources 
brighter than flux 5* is Nsrc{> S), while the survey flux limit 
is 5*0. The survey is sensitive to systems in which the in- 
trinsic flux S and the lensing magniflcation ^ combine such 
that fj,S > So- Thus, the total number of lenses found in the 
survey has the form of an integral over the source plane, 



iMitchefl et~aLll2005l : iHuterer et al.ll2005h : 



Mens OC / iVsrc(>5'o//i) d/3. 



(35) 



fTn = / 

J n 



Nsrc{>So/fi) 
Nsrc{>So) 



d/3 



(36) 



The subscript n indicates that we can compute the biased 
cross-section for lenses with n images by restricting the in- 
tegral to the appropriate region in the source plane (cf. 
Fig. I10|l . (Note that we do not distinguish between the sym- 
bols for unbiased and biased cross-sections, because we have 
tried to be explicit in the subsequent text about which we 
are using.) The cross-section can be combined with the num- 
ber density of sources to obtain the observed number of lens 
systems. 

While magniflcation bias depends explicitly on the flux 
limit of a survey, it depends implicitly on the resolution as 
well. For a survey that detects lenses but does not resolve 
the images (assuming they are resolved in follow-up observa- 
tions) , the initial detection depends on the combined flux of 
all images, so it is the total magniflcation that enters equa- 
tion (|36|) . For a survey that can resolve the images at the 
outset, the discovery of lensing depends on the flux of the 
second-brightest image so it is the second-ranked magni- 
flcation that is used in equation (|36|l . To keep our analysis 
as general as possible, we present results for both cases. 

To compute biased cross-sections we must specify the 
source counts, Nsrc{> S), or equivalently the differential 
counts, dNsrc/dS, near the survey flux limit. We consider 
optical quasar surveys to determine a reasonable model. 
For optical surveys, dN/dS is typically modelled as a bro- 
ken power-law, so we consider the two extremes. Based 
on the results for S PSS and 2QZ in the g- and i-bands 
jCroom et all |2004| : [Richards fc et al.1 12OO6I I , we flnd that 
dlog A^src/dmag = a where a ~ 1 at the bright end and 
~ 0.1 at the faint end. This corresponds to 

^oc5 , (37) 

where v — 2.5a + 1. Thus, our two values of 1/ are 3.5 and 
1.25. (A faint-end slope of = 1.25 is consistent with the 
spectroscopic survey of faint quasars by Jiang & et al. 200ti.) 
To account for the break in the power-law, we implement a 
functional form commonly used for quasar luminosity func- 
tions. 



$(L) (X 



(38) 



This number is divided by the total number of sources to ob- 
tain the lensing probability, lead ing to the following deflni- 
tion for the biased cross-section (jKeeton fc Zabludofill20o4 



where we use vi and 1^2 of 3.5 and 1.25. For simplicity, we 
consider limiting luminosities rather than limiting fluxes, 
but we examine a broad range of values: 0.04L, , 0.4L« , and 
4L, (chosen deliberately to sample the full range of effective 
slopes with even spacing in magnitude). These correspond 
to effective slopes VcS in equation (|37|l of 1.25, 1.5, and 3.4, 
although we caution that the luminosity function is far from 
a pure power-law (especially near the intermediate limiting 
luminosity), so magniflcation bias depends not just on ly^B 
but on the full shape of the luminosity function. 

A 2-image system in which one image is fainter than the flux 
Hmit will not be identified as a lens. A system with more than 
two images will be identified as a lens if at least two images are 
above the flux limit. (If some images are below the flux limit, the 
true number of images may not be known at flrst, but the system 
will still be identifled as a lens.) 
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Figure 10. Sample mock lenses for our very oblate higher mass, group scale model viewed edge-on. The left panel shows the caustics 
in the source plane, along with random source positions selected with magnification weighting (see Section |5. 4^ . The right panel shows 
the critical curves in the lens plane, along with the images of the random sources. The horizontal axis is aligned with the major axis of 
the surface mass density of the dark matter component (which in turn is aligned with the stellar component in this axisymmetric case). 
The point type denotes the number of images, as indicated in the legend. The left panel clearly shows the division of the source plane 
into 2, 3, and 4-image regions bounded by the caustics. The right panel shows that the lens plane is also divided into distinct regions, 
although they are not bounded by easy-to-find curves (cf. Fig. 1 of lPinch et alll2002l) . This figure includes only 2000 sources for clarity, 
compared with 5000 used in our full analysis. 



In summary, we track 7 values of cross-section: the un- 
biased one, and 6 biased ones (= 2 weighting schemes x 3 
limiting magnitudes). The GRAVLENS routine mockLF (new 
in version 1.99k) can be used to define a double power-law 
source luminosity function that mockcsec uses to compute 
biased cross-sections. 



5.5 Image separation distributions 

In addition to the lensing cross-sections, we also consider 
the distribution of image separations A9. Here, we define 
the image separation for a given lensing system as the max- 
imum distance between any two images in the system. This 
definition is simple and well-defined even for image sys- 
tems with arbit rary numbers of images, and it has been 
used before (e.g..lliiada fc et al.ll2003l : lOguri fc Keetonl2004 



iHuterer et al.ll2005l ). 

In principle, especially for systems with > 2 images, 
there may be more sophisticated ways to characterise im- 
age separations that incorpor ate informa t ion ab out t he full 
config uration of images. While iKochane^ l|2002b[ ) and lOguril 
l|2007f ) have taken first steps in this direction, there has been 
no systematic exploration of various image separation statis- 
tics. We defer consideration of these issues to future work. 



radius to be Rain ~ 1" and we consider a simple power-law 
source luminosity function with v = 1.25, as appropriate 
for the faint end of our fiducial broken power-law model. 
(The biased cross-section cannot be computed analytically 
for the full broken power-law luminosity function.) The ratio 
of the numerical to analytic cross-section is 0.988 ±0.003 for 
the unbiased cross-section, as well as for the biased cross- 
section using the total magnification. (The statistical er- 
ror bar is from bootstrap resampling.) Thus, the numerical 
cross-sections are accurate to about 1 per cent. 

To test our Fourier methods for handling map-based 
lens models, we present an example using our higher mass 
scale deprojected Sersic model. We compute the cross- 
sections with an analytic circular lens model, now using our 
fiducial broken power-law source luminosity function and 
tracking 7 cross-sections (one raw and six biased). We then 
create a /t-map that is 30" on a side and has 800^ pixels 
and recompute the cross-sections using a map-based lens 
model. The ratio of the cross-section from the map-based 
model to that from the analytic model is between 0.960 and 
1.041 for all 7 cross-sections. The statistical uncertainties 
depend on the bias mode, ranging from 0.001 (for the un- 
biased cross-section) to 0.023 (for the biased cross-section 



5.6 Validation 

To validate our Monte Carlo methods for computing lensing 
cross-sections, we present a simple test case with a singu- 
lar isothermal sphere lens, for which cross-sections can be 
computed analytically. For concreteness we set the Einstein 



^ The box size is the fiducial value determined from the box 
size tests discussed in Section 15.81 and the resolution is twice 
the fiducial resolution discussed in Section 15.71 While we clearly 
show in that subsection that the fiducial resolution is sufficient 
for cusped models, we have found that twice better resolution is 
necessary to resolve the very steep inner parts of the higher mass, 
deprojected Sersic model (stellar component). 
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using the second-brightest magnification and a limiting lu- 
minosity of 4L,). The ratios become even closer to unity as 
the resolution increases, but we defer further discussion of 
the resolution to Section [5]7] The main conclusions we draw 
from this test are that the Fourier methods are valid and 
that map-based lens models yield accurate cross-sections. 

5.7 Resolution 

As noted in Section [5.11 we need to ensure that the lensing 
properties of our model galaxies are not affected by the res- 
olution of the grid on which we solve the Poisson equation. 
Here we present convergence tests to determine acceptable 
pixel sizes. Before designing these tests, we must determine 
precisely how the pixel size can cause problems with res- 
olution. There are two separate issues: (i) In cases where 
the Einstein radius is quite small, the finite pixel size may 
cause GRAVLENS to calculate a polygonal rather than ellipti- 
cal inner critical curve (in the lens plane), and consequently 
a polygonal outer source plane caustic as well, (ii) In cases 
where the inner slope of the density profile is very steep, a 
large pixel size may mean that its slope, and therefore the 
location of the inner critical curve, are poorly determined. 
We consider the two cases separately below. 

5.7.1 Obtaining a smooth inner critical curve 

For these tests, we use the mass models with 7dm = 0.5, since 
they have the smallest Einstein radii and cross-sections, and 
therefore yield the most stringent resolution requirements 
from a standpoint of having a smooth inner critical curve. 

We begin with a pixel size of 0.025" and 0.075" for the 
galaxy and group scale models, respectively (with 10" and 
30" boxes; see Section[5]8}. We vary the pixel size to a factor 
of 2 smaller, and a factor of 2 and 4 larger in the hopes (a) of 
seeing that the lensing results have converged by the smallest 
pixel sizes, and (b) of determining the largest pixel sizes that 
yield reliable results. We do these tests using a limited set of 
models: a spherical model, followed by the very oblate, very 
prolate and triaxial models seen from the most extreme 2, 
2, and 3 viewing directions along the principal axes. 

The figures of merit for the convergence tests are the 
lensing cross-sections (cr2, crs, 0-4 and atot separately), al- 
though we also examine the critical curves and caustics. 
We are interested in cross-sections both with and without 
magnification bias (according to the various magnification 
schemes discussed in Section 15. 4[) . For the purpose of this 
test, we require convergence of the cross-sections at the 5 
per cent level. 

The effects of resolution are seen qualitatively in the 
caustics, which are shown in Fig. [TT] for three of the eight 
lens configurations. (Note that the mass density distribu- 
tions are always elongated along the x' axis, but as usual 
the caustics are elongated in the perpendicular direction.) 
The quantitative results of the resolution tests can be sum- 
marised as follows: 

• For the lower mass scale, the most stringent require- 
ments on the resolution come from the very oblate, edge- 
on model, for which we require 0.025" pixels. Even for this 
model, the resolution seems to have converged, i.e., the dif- 
ference between 0.025" and 0.0125" pixels is much less than 
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Figure 11. Source plane caustics for the lower (left) and higher 
(right) mass models with three of the eight extreme lens configu- 
rations used for the '^^^ = 0.5 resolution tests: spherical, oblate 
(edge-on), and prolate (side-on) from top to bottom. The pixel 
scales are as labelled on the top plot for each mass scale. The axes 
are in the plane of the sky (in units of arcsec) with the horizontal 
axis aligned with the major axis of the surface mass density of 
the dark matter component. 



our 5 per cent convergence criterion. For the very prolate, 
side-on model, 0.025" pixels are formally insufficient for the 
4-image cross-section. However, the 4-image cross-section is 
only ~ 5 per cent of the total cross-section, and with such 
a small value the statistical uncertainties are large, so the 
failure is actually marginal and not a significant concern. 

• For the higher mass scale, the most stringent resolu- 
tion test also comes from the very oblate, edge-on case. 
The second-best resolution (0.075" pixels) marginally sat- 
isfies our convergence criterion for the 4-image case, and 
marginally fails for the 3-image case. However, this marginal 
failure is unlikely to be too crucial, since the 3-image cross- 
section is such a small fraction of the total cross-section. In 
all other cases, the second-best resolution (0.075" pixels) is 
clearly sufficient at the 5 per cent level. 

• In general, as we have predicted, the radial caustic is 
much more strongly affected by pixelisation effects than the 
tangential caustic, since the radial caustic maps to the inner 
critical curve in the lens plane. The radial caustic tends to 
become more polygonal rather than elliptical as the resolu- 
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tion is degraded. In Fig. 1111 it is clear that the caustics are 
nearly identical for both the best and second-best resolu- 
tions, and that the two poorer resolutions are insufficient. 

In conclusion, our second-best resolutions for both mass 
scales , 0.025" (~ 0.11 kpc) and 0.075" (~ 0.33 kpc), appear 
to be safe choices from the standpoint of obtaining smooth 
critical curves and caustics, and the results have converged 
in comparison with higher resolution simulations at the 5 
per cent level. In dimensionless units, these resolutions cor- 
respond to ~ 25 pixels lengths per -Rein. 

5.7.2 Resolving a steep inner slope 

The other resolution issue, the ability to resolve the slope of 
the density profile in the inner regions, is most easily tested 
using our steepest cusped model, 7dm ~ 1.5. In this case, we 
repeat the previous tests but only for the spherical case. 

Our results suggest that while this efi'ect may be sig- 
nificant for the poorest resolution we consider, the results 
have converged at better than the 5 per cent level for our 
fiducial resolution. For example, in the previous subsection 
(spherical 7dm = 0.5 case), the fiducial resolution gave cross- 
sections that were 1.4 and 2.7 per cent lower than the best 
resolution (lower and higher mass models, respectively); the 
poorest resolution led to a highly significant 20 and 35 per 
cent reduction of the cross-section. For 7dm ~ 1.5, those 
numbers are 1.1 and 0.9 per cent reduction in atot for the 
fiducial resolution, or 15 and 11 per cent reduction for the 
poorest resolution. In all cases, these results are for the unbi- 
ased cross-section, for which it is clear that the fiducial reso- 
lution is sufficient to solve both possible resolution problems. 
For the biased cross-section, the resolution requirements ap- 
pear to be slightly less stringent. 

While we do not use a map-based approach for the 
deprojected Sersic profile simulations, similar convergence 
tests suggest that such an approach would lead to more 
stringent resolution requirements (at least a factor of two 
better), due to the necessity of resolving a very steep inner 
slope in the innermost regions. The reason these models may 
be particularly problematic is that they do not converge to 
a single inner slope, unlike the cusped models. 

5.8 Box size 

We also need to check that the size of the box in which 
we solve the Poisson equation does not aff'ect our lensing 
results. We have again run convergence tests, but now using 
mass models with 7dm ~ 1.5 since these have the largest 
Einstein radii and hence the greatest sensitivity to the box 
size. We include non-spherical models since the elongation 
of the mass distribution may create extra demands on the 
size of the box in the direction of the major axis. As in 
the resolution tests, we use only the extreme non-spherical 
configurations, as a worst-case scenario. 

To begin, we use 10" x 10" and 30" x 30" boxes for 
the galaxy and group scale models, respectively (with our 
fiducial 0.025" and 0.075" pixels; see Section [53. We then 
examine box sizes of 5", 10", and 15" for galaxy scale mod- 
els, and 15", 30", and 45" for group scale models. Fig. 1121 
shows the efi^ects of the box size on the caustics, for three 
of the eight extreme lens configurations. A summary of the 
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Figure 12. Source plane eausties for the lower (left) and higher 
(right) mass models, with the three of the eight extreme lens con- 
figurations used for the box size tests: spherical, edge-on oblate, 
and side-on prolate (top to bottom) . Box sizes are labelled on the 
top plot for each mass scale. The meaning of the axes is the same 
as in Fig. [11] 



quantitative results for all eight lens configurations are as 
follows: 

• For the lower mass scale, the default box size of 10" 
is sufficient. There is one marginal case — the 3-image cross- 
section for the very oblate, edge-on model — but this cross- 
section is a small fraction (< 1 per cent) of the total cross- 
section for that model, and the statistical errors in deter- 
mining it are large. The smallest (5") box typically yields 
unbiased cross-sections that pass our 5 per cent convergence 
criterion, but the biased cross-sections can fail dramatically 
(up to 25 per cent reduction of atot relative to the largest 
box), which eliminates the possibility of using such a small 
box. 

• For the higher mass scale, the default box size of 30" 
is typically sufficient to achieve 5 per cent accuracy com- 
pared to the larger box size, though occasionally the results 
are marginal, particularly for highly biased cross-sections. 
Visually, the caustics appear to be nearly identical for the 
30" (fiducial) and 45" (largest) boxes (Fig. [12}, with rather 
bizarre and very significant distortion for the 15" box for 
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some viewing directions, for which the cross-sections are also 
significantly reduced. 

• Some caustics appear smoother than for the corre- 
sponding cases in the resolution tests in Section TS.?! The dif- 
ference comes from the fact that here we use 7dm = 1.5 mod- 
els, which have Einstein radii that are significantly larger 
than for the 7dm = 0.5 models used for the resolution tests. 

• The effect of insufficient box size is to shrink both the 
radial and tangential caustics, and to alter the shape of the 
tangential caustic. 

The caustics and cross-sections together imply that even 
with the largest possible Einstein radius at each mass 
scale, the medium-sized boxes, 10" (~ 43.6 kpc) and 30" 
(~ 131 kpc), are sufficient. The box size requirement in di- 
mensionless units is ~ 16i?cin. 

5.9 Sampling viewing angles 

We are interested in the lensing properties of galaxies not 
only for specific viewing angles, but also when averaged over 
all viewing angles. We must choose how to sample the pos- 
sible viewing directions in order to obtain accurate averages 
with as little computational effort as possible. We are guided 
by two principles. First, we expect the lensing quantities 
to vary with viewing direction in some smooth and, in the 
oblate and prolate cases, monotonic way. Second, this vari- 
ation is more directly related to the shape of the projected 
surface mass density than to the viewing angles per se. 

5.9.1 Sampling in projected axis ratios 

Instead of directly sampling angles relative to the intrin- 
sic axes -& and ip, we sample linearly in the ratios a' /a and 
fe'/a, i.e., the projected major and minor semi-cixis length 
a' and b' as given in equation ((S}, normalised by the in- 
trinsic length scale a. Since this choice does not correspond 
to an equal-area sampling on the viewing sphere (i.e., uni- 
form in — 1 ^ cosi9 ^ 1 and Q ^ tp ^ 27r), we need to 
include a weighting factor to recover accurate averages. The 
orientation-averaged value of some quantity / (which in our 
case will be a lensing cross-section) can be written as 

(/) = - r^'/(^,^) sin^d^dy., (39) 

1" Jo Jo 




f {d, if) sin-d J{a ,b') Aa Ab' , (40) 



where the expression for the Jacobian J(a', b') is given in 
AppendixlA] as is the inversion of equation ((5)1 to find (i9, if) 
in terms of (a', 6'). In the oblate and prolate cases, one of 
the two integrals is irrelevant since the lensing depends on 
only one of the angles, but in the triaxial case both integrals 
are important. 

By construction, the individual mass components we 
use have axis ratios that are constant as a function of ra- 
dius. However, since the dark matter and stellar components 
have different intrinsic shapes (Tabled, the combined, total 
density has intrinsic and projected axis ratios that vary with 
radius (see also Fig. |9]). Fortunately, the precise choice of 
axis ratios to use for the sampling is not terribly important, 
as long as it is treated properly according to equation (|40p . 



We choose to sample in the axis ratios of the dark matter 
component, i.e., in what follows a' /a and b' /a always refer 
to the dark matter component. 

There are two important issues related to the sampling: 

(i) We must first test our intuition that cross-sections 
vary smoothly with a' /a and b' /a. We can do this by using a 
large number A'^s of samplings (we use Ns = 47 for the oblate 
case, Ns = 31 for the prolate case, and Ns = 23 x 15 = 345 
for the triaxial case). 

(ii) We must determine how few samplings can be used to 
obtain orientation-averaged cross-sections that are accurate 
at the 5 per cent level. We can do this by undersampling the 
original number of samplings and seeing how few samplings 
we can use while still recovering the average cross-sections 
to within 5 per cent. 

For the viewing angle sampling tests, we use cusped 
density profiles with dark matter slope 7dm = 1.5 and the 
extreme oblate and prolate shapes, since these will have the 
strongest variation of cross-section with inclination. We also 
consider the triaxial case, for which we test the sampling in 
both a' /a and fo'/a. In addition to the total strong lensing 
cross-section, we require accurate recovery of the 2, 3, and 
4-image cross-sections separately. 

5.9.2 Variation of cross-section with projected axis ratios 

We show the variation of the cross-sections with viewing 
angle in Fig. [13] for both mass models, using the very oblate 
and very prolate shapes. The symmetry of the oblate case 
means we only need to vary the polar viewing angle i?, which 
we accomplish by taking equal steps in a' /a. Similarly, for 
the prolate case we only need to vary the azimuthal viewing 
angle ip, again taking equal steps in a' /a. For this plots, we 
use our fiducial number of samplings (47 for the oblate case, 
and 23 for the prolate case) to generate the cross-section 
curves, and we show the undersampled by four case with 
crosses. We defer discussion of the physical intuition behind 
the viewing angle dependence to Paper II, and focus here 
on the technical issues of smoothness and sampling. It is 
clear that, as expected, the variation of the individual and 
total cross-sections is a smooth function of the projected 
axis ratios for prolate and oblate dark matter plus stellar 
density shapes. 

Note that the 3-image cross-section should be identi- 
cally zero for any lens configuration for which the tangential 
caustic lies entirely within the radial caustic. This means 
that CT3 should be zero for our oblate models when fe'/a is 
above some threshold, and for our prolate models through- 
out the full range of a' /a. Our numerical 3-image cross- 
sections are not strictly zero at a very few values of the 
projected axis ratio due to extremely rare errors in finding 
and classifying images (see Section [5]2]). However, the spuri- 
ous 3-image cross-sections are well below 0.1 per cent of the 
total cross-sections, so they are a negligible source of error. 

In Fig. 1141 we show the variation of the unbiased cross- 
sections with a' /a and fe'/a for the lower mass triaxial model. 
(The results are very similar for the higher mass triaxial 
model, except for an overall difference in amplitude.) The 
changes for the biased cross-sections are described qualita- 
tively. We again defer discussion of the physical intuition be- 
hind the orientation-dependence of the cross-section to Pa- 
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Lower mass scale, oblate (left) and prolate (right) 



Higher mass scale, oblate (left) and prolate (right) 
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Figure 13. Lensing cross-sections for the lower (left plot) and higlier (riglit plot) mass scale, for oblate and prolate models, for 2, 
3, 4-image systems separately as labelled, and the total cross-section at the bottom. The lines are as follows: black solid = unbiased, 
red dotted and blue dashed = weighted by total and second-brightest magnification respectively for 0.04L, limiting magnitude, green 
long-dashed and yellow dot-dashed = the same for 0.4L» limiting magnitude, and magenta dot-long dashed and cyan short-long dashed 
= the same for 4L* limiting magnitude. The lines show the curves resulting from our fiducial sampling in o! and 6'; crosses show our 
eventual choice of sampling as described in Section 15.9.31 The 3-image cross-section should be identically zero for fe'/a or a'/a above 
some threshold; the occasional spikes at small but finite values indicate small numerical errors due to misclassification of images. 



per II. Here we merely observe that for the triaxial case, the 
unbiased cross-sections a"2, o"4, and crtot are smooth functions 
of a'/a and h' ja, without multiple local minima/maxima or 
other features that would require dense sampling. The bi- 
ased cross-sections also vary smoothly, regardless of which 
magnification bias mode we use. One difference is that the 
contours of the 2-image biased cross-section tend to lie at 
nearly constant a'/a, not as tilted as in the unbiased case 
(see Paper II). 

5.9.3 Minimum number of samplings 

Next we consider how few samplings we can get away with 
and still recover accurate results for cross-sections averaged 
over viewing directions. To do so, we start with our fidu- 
cial number of samplings, and compare the total average 
cross-sections (properly weighting each sampling using the 
Jacobian) when undersampling by a factor of 2, 4, and 8. 
In the triaxial case, we always consider undersampling by 
the same factor in both a'/a and b'/a. Generally, given Na 
samplings in a'/a or b'/a, we have {1 + Ns) /u viewing angles 
when undersampling by a factor of = 2" for « J5 1. 

We find that for all mass scales, magnification bias 
modes, and numbers of images, undersampling by a factor 
of 4 in both a'/a and b'/a allows us to recover orientation- 
averaged cross-sections to within 5 per cent accuracy when 
compared with the full sampling. In general, 3-image sys- 
tems are the most difficult to sample properly due to their 
strong variation with inclination (as in Fig. I13p . This level 



of undersampling corresponds to only 12 samplings in b'/a 
in the oblate case, 8 samplings in a'/a in the prolate case, 
and 6 (in a'/a) x 4 (in b'/a) = 24 samplings total in the 
triaxial case. For reference, in Fig.ll3land ll4l we indicate our 
chosen level of undersampling with crosses. 



6 CONCLUSIONS 

We have presented a fiexible simulation pipeline for coher- 
ent investigations of selection and modelling biases in strong 
lensing surveys. We have focused on point-source lensing 
by two-component galaxy models meant to emulate realistic 
early- type, central galaxies at two different mass scales: a 
lower, ~ 2L, galaxy mass scale and a higher, ~ 7Lt group 
mass scale (with L^, in the r-band). Below is a list of our 
main conclusions regarding the construction of the simula- 
tion pipeline for lensing by realistic galaxy models: 

• We include both cusped and deprojected Sersic density 
profiles, with observationally-motivated choices of masses 
and scale lengths, and a range in density profile parame- 
ters, separately for the stellar and dark matter component. 
We use seven different models for the galaxy shapes with the 
stellar component rounder than the dark matter component, 
but with the axes intrinsically aligned. [Sections 14 . 1144 . 5] 

• When we change the density profile parameters away 
from the adopted fiducial values, we preserve the total and 
virial mass of the stellar and dark matter components (re- 
spectively) by changing the density amplitude. We keep the 
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Figure 14. Base-10 logarithm of the unbiased lensing cross- 
sections for the lower mass, triaxial shape model, as a function of 
a' /a and fe'/a. Note that the contour scales are different on each 
plot, and (73 is not shown because it is negligible. The crosses in- 
dicate our final choice of sampling, as described in Section 15.9.31 



scale radius fixed to comply with observed size-luminosity 
relations for early-type galaxies and concentration-mass re- 
lations for dark matter halos. Similarly, for the non-spherical 
shapes we change the scale length to preserve the mass but 
keep the scale radius (and hence the concentration) the same 
as for the fiducial spherical case. [Sections 14. 3. 21 H. 4. 21 14. 5."^ 

• We choose fiducial lens and source redshifts of 0.3 and 
2.0, resulting in Einstein radii Rdn at about one effective 
radius Re- [Section 13.2) 

• We require multiple ways of handling magnification bias 
to sample different parts of the observed quasar luminosity 
function and to allow for application to different survey lim- 
itations in image resolution and flux. [Section |5. 4) 

• When the surface mass density is known analytically 
(as for our deprojected Sersic models), we compute the cor- 
responding lensing deflection and magnification analytically 
(for circular symmetry) or with numerical integrals (for el- 
liptical symmetry). When the surface mass density cannot 
be computed analytically, as for the cusped models, we con- 
struct a surface density map and use Fourier methods to 
compute the lensing defiection and magnification. We val- 
idate and use Monte Carlo methods to calculate lensing 
statistics, including unbiased and biased cross-sections and 
image separations. [Sections 15. 1H5. 6] 



• For our map-based lensing calculations, tests for conver- 
gence at the 5 per cent level suggest that the required map 
resolution is ~ 25 pixel lengths per -Rein, corresponding to 
0.025" (~ 0.11 kpc) and 0.075" (~ 0.33 kpc) for the lower 
and higher mass scale, respectively. The resolution needs to 
be fine enough to simultaneously resolve the steepness of the 
inner slope of the density and recover the smoothness of the 
shape of the inner critical curve. Similarly, convergence tests 
suggest that the required box size is ~ 16 -Rein, correspond- 
ing to 10" (~ 43.6 kpc) and 30" (~ 131 kpc) for the lower 
and higher mass scale, respectively. [Sections 15.71 and I5.8j 

• Lensing cross-sections for fairly general triaxial models 
vary in a smooth and monotonic way with viewing direction, 
which can be efficiently sampled through projected axis ra- 
tios (rather than viewing angles) of the surface density. Sur- 
prisingly few samplings axe necessary: 10, 8, and 24 viewing 
directions can effectively determine the orientation-averaged 
cross-sections for oblate, prolate, and triaxial models (re- 
spectively) with 5 per cent precision. [Section |5. 9) 

We have implemented the construction of the realistic galaxy 
models and efficient computation of the corresponding sur- 
face mass density maps in idl (see also Section 14. 6p . All 
subsequent strong lensing calculations are done with an up- 
dated and extended version of gravlens. We use scripts to 
coordinate and analyse the extensive data flow to and from 
these two parts, and to make it into a coherent simulation 
pipeline. 

Our investigations have revealed a number of useful, 
general points related to the interpretation of observational 
strong lensing results: 

• In the vicinity of the Einstein radius, the intrinsic and 
projected logarithmic slopes of the total density are close to 
the values 7 = 2 and 7' = 1 of an isothermal profile, for 
all of our galaxy models (see Figs. [3] and O and the text of 
Sec. l4.3.3|l . Consequently, a measurement of the total density 
profile near the Einstein radius cannot (alone) be used to 
determine the dark matter inner slope for a two-component 
model. This is a non-trivial conclusion given the variety of 
intrinsic non-isothermal density profiles used for both the 
stellar and dark matter component, including profiles with 
a true central cusp and those that do not asymptote to a 
particular inner slope at any scale. 

• We also find that the lensing defiection curves are nearly 
flat around (and even beyond) -Roin, similar to a flat deflec- 
tion curve for an isothermal model (see Fig. (Tjl , particularly 
for the lower-mass galaxy scale. This again implies that con- 
straints from strong lensing cannot be extrapolated to radii 
much smaller or larger than the Einstein radius. We empha- 
size that the parameters for the galaxy models were not a 
priori chosen to mimic this effect, b ut that it truly see ms to 
bean "isothermal conspiracy" fe.g.. lRusin et all(2003l ). 

• Because of our assumption of alignment between the 
intrinsic axes of the stellar and dark matter components, 
only the two shape models with a triaxial dark matter com- 
ponent and a rounder (triaxial or oblate) stellar component 
can have a significant misalignment between their projected 
axes. Even then, significant misalignment occurs only for a 
limited number of viewing directions, several of which lead to 
a relatively round projected stellar component such that the 
misalignment is not very important in practice (the position 
angle of a round stellar component is on the sky, after all. 
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difficult to establish). The small misalignments inferred for 
relatively isolated lens systems, as well as detailed dynam- 
ical modelling of early-type galaxies, support near intrinsic 
alignment between stars and dark matter (see Section|4A3|. 

In Paper II, we use the flexibility of the pipeline to 
study selection biases related to the galaxy mass, shape, 
orientation, and various parameters of the dark matter and 
stellar profiles. In subsequent work, we will investigate mod- 
elling biases by analysing the mock lens systems produced 
by the pipeline using the lens modelling tools within the 
GRAVLENS package. We have therefore created a lensing sim- 
ulation pipeline that will be crucial for deriving inferences 
about galaxy density profiles and shapes, Hq, and other pa- 
rameters from the thousands of lenses that will be discov- 
ered in the coming years in large photometric surveys such 
as Pan-STARRS, LSST, SNAP, and SKA. 

Furthermore , following a similar approach as in 
Ivan de Ven et ahl (2008a), we plan to extend the pipeline 
to produce projected kinematics of our galaxy models that 
mimic the two-dimensional kinematic observations of early- 
type galaxies. We will then investigate how well we can ex- 
pect to recover the intrinsic density profile and shape of 
early-type galaxies, based solely on dynamical models fit- 
ted to these simulated kinematics, or on a combination of 
kinematics and s trong lensing (using techniques similar to 
Ivan de Ven et ahl 2008b). This work will be valuable for un- 
derstanding the modelling biases in analyses of the kine- 
matic data that are becoming available for hundreds of 
galaxies at increasing redshifts, and for understanding both 
modelling and selection biases in joint lensing+kinematics 
studies. 
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APPENDIX A: SAMPLING OF VIEWING 
DIRECTION VIA PROJECTED AXIS RATIOS 

The smooth variation in lensing cross-section with viewing 
direction is more directly related to the (elliptic) shape of 
the surface mass density than the viewing angles. Hence, 
instead of sampling directly ^ and ip, we sample linearly in 
the projected axis ratios a' and b' given in equation 
Inversion of the latter relation yields 

cos'^ - (a ~c){b -c) 

tan (z; = -^^ 4^ '-. (A2 

The Jacobian J{a',b') of the transformation between {"ffjip) 
and (a', b') is given by 

. ojf'u'. a'b'(a'^-b'''){b''~c^) 
smt>J(a ,b)- _ ^^^^^^ _ ^,2^^^,, _ ^^^^^,2 _ • 

(A3) 

In the axisymmetric limit one of the viewing angles becomes 
redundant. In the oblate case (a = & > c), we set, without 
loss of generality, ip = n /2 so that a' = a, and 

cos^-ff^— -, siTn9H79= (A4 

a^-c^ 7(a2 _c2)(6'2 _c2) 

Similarly, in the prolate case (a > fe = c) we take = n/2 
so that b' = c, and 

2 a'^ — , o,' da' , , 

tan (p = ^- ^, d<z) = — A5) 

^ a2_„,2' ^ ^(a2_a'2)(a'2-c2) ^ ^ 
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